Interpretation and design of hydraulic fracturing treatments

ABSTRACT

Solutions for the propagation of a hydraulic fracture in a permeable elastic rock and driven by injection of a Newtonian fluid. Through scaling, the dependence of the solution on the problem parameters is reduced to a small number of dimensionless parameters.

RELATED APPLICATION

This application is a continuation of U.S. patent application Ser. No.10/356,373, filed Jan. 31, 2003, now U.S. Pat. No. 7,111,681 whichclaims the benefit of priority under 35 U.S.C. 119(e) to U.S.Provisional Patent Application Ser. No. 60/353,413, filed Feb. 1, 2002,which applications are incorporated herein by reference.

FIELD

The present invention relates generally to fluid flow, and morespecifically to fluid flow in hydraulic fracturing operations.

BACKGROUND

A particular class of fractures in the Earth develops as a result ofinternal pressurization by a viscous fluid. These fractures are eitherman-made hydraulic fractures created by injecting a viscous fluid from aborehole, or natural fractures such as kilometers-long volcanic dikesdriven by magma coming from the upper mantle beneath the Earth's crust.Man-made hydraulic fracturing “treatments” have been performed for manydecades, and for many purposes, including the recovery of oil and gasfrom underground hydrocarbon reservoirs.

Despite the decades-long practice of hydraulic fracturing, manyquestions remain with respect to the dynamics of the process. Questionssuch as: how is the fracture evolving in shape and size; how is thefracturing pressure varying with time; what is the process dependence onthe properties of the rock, on the in situ stresses, on the propertiesof both the fracturing fluid and the pore fluid, and on the boundaryconditions? Some of the difficulties of answering these questionsoriginate from the non-linear nature of the equation governing the flowof fluid in the fracture, the non-local character of the elasticresponse of the fracture, and the time-dependence of the equationgoverning the exchange of fluid between the fracture and the rock.Non-locality, non-linearity, and history-dependence conspire to yield acomplex solution structure that involves coupled processes at multiplesmall scales near the tip of the fracture.

Early modeling efforts focused on analytical solutions for fluid-drivenfractures of simple geometry, either straight in-plane strain orpenny-shaped. They were mainly motivated by the problem of designinghydraulic fracturing treatments. These solutions were typicallyconstructed, however, with strong ad hoc assumptions not clearlysupported by relevant physical arguments. In recent years, thelimitations of these solutions have shifted the focus of research in thepetroleum industry towards the development of numerical algorithms tomodel the three-dimensional propagation of hydraulic fractures inlayered strata characterized by different mechanical properties and/orin-situ stresses. Devising a method that can robustly and accuratelysolve the set of coupled non-linear history-dependentintegro-differential equations governing this problem will advance theability to predict and interactively control the dynamic behavior ofhydraulic fracture propagation.

For the reasons stated above, and for other reasons stated below whichwill become apparent to those skilled in the art upon reading andunderstanding the present specification, there is a need in the art foralternate methods for modeling various behaviors of hydraulic fracturingoperations.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a view of a radial fluid-driven fracture with anexaggerated aperture;

FIG. 2 shows a tip of a fluid-driven fracture with lag;

FIG. 3 shows a rectangular parametric space;

FIG. 4 shows a pyramid-shaped parametric space;

FIG. 5 shows a triangular parametric space;

FIG. 6 shows a semi-infinite fluid-driven crack propagating in elastic,permeable rock;

FIG. 7 shows another triangular parametric space;

FIG. 8 shows a plane strain hydraulic fracture;

FIG. 9 shows another rectangular parametric space;

FIG. 10 shows a triangular parametric space with two trajectories;

FIG. 11 shows a graph illustrating the dependence of a dimensionlessfracture radius on a dimensionless toughness; and

FIG. 12 shows another triangular parametric space with two trajectories.

DESCRIPTION OF EMBODIMENTS

In the following detailed description, reference is made to theaccompanying drawings that show, by way of illustration, specificembodiments in which the invention may be practiced. These embodimentsare described in sufficient detail to enable those skilled in the art topractice the invention. It is to be understood that the variousembodiments of the invention, although different, are not necessarilymutually exclusive. For example, a particular feature, structure, orcharacteristic described herein in connection with one embodiment may beimplemented within other embodiments without departing from the spiritand scope of the invention. In addition, it is to be understood that thelocation or arrangement of individual elements within each disclosedembodiment may be modified without departing from the spirit and scopeof the invention. The following detailed description is, therefore, notto be taken in a limiting sense, and the scope of the present inventionis defined only by the appended claims, appropriately interpreted, alongwith the full range of equivalents to which the claims are entitled. Inthe drawings, like numerals refer to the same or similar functionalitythroughout the several views.

The processes associated with hydraulic fracturing include injecting aviscous fluid into a well under high pressure to initiate and propagatea fracture. The design of a treatment relies on the ability to predictthe opening and the size of the fracture as well as the pressure of thefracturing fluid, as a function of the properties of the rock and thefluid. However, in view of the great uncertainty in the in-situconditions, it is helpful to identify key dimensionless parameters andto understand the dependence of the hydraulic fracturing process onthese parameters. In that respect, the availability of solutions foridealized situations can be very valuable. For example, idealizedsituations such as penny-shaped (or “radial”) fluid-driven fractures andplane strain (often referred to as “KGD,” an acronym from the names ofresearchers) fluid-driven fractures offer promise. Furthermore, the twotypes of simple geometries (radial and planar) are fundamentally relatedto the two basic types of boundary conditions corresponding to the fluid“point”-source and the fluid “line”-source, respectively.

Various embodiments of the present invention create opportunities forsignificant improvement in the design of hydraulic fracturing treatmentsin petroleum industry. For example, numerical algorithms used forsimulation of actual hydraulic fracturing treatments in varying stressenvironment in inhomogeneous rock mass, can be significantly improved byembedding the correct evolving structure of the tip solution asdescribed herein. Also for example, various solutions of a radialfracture in homogeneous rock and constant in-situ stress presentnon-trivial benchmark problems for the numerical codes for realistichydraulic fractures in layered rocks and changing stress environment.Also, mapping of the solution in a reduced dimensionless parametricspace opens an opportunity for a rigorous solution of an inverse problemof identification of the parameters which characterize the reservoirrock and the in-situ state of stress from the data collected duringhydraulic fracturing treatment.

Various applications of man-made hydraulic fractures includesequestration of CO₂ in deep geological layers, stimulation ofgeothermal reservoirs and hydrocarbon reservoirs, cuttings reinjection,preconditioning of a rock mass in mining operations, progressive closureof a mine roof, and determination of in-situ stresses at great depth.Injection of fluid under pressure into fracture systems at depth canalso be used to trigger earthquakes, and holds promise as a technique tocontrol energy release along active fault systems.

Mathematical models of hydraulic fractures propagating in permeablerocks should account for the primary physical mechanisms involved,namely, deformation of the rock, fracturing or creation of new surfacesin the rock, flow of viscous fluid in the fracture, and leak-off of thefracturing fluid into the permeable rock. The parameters quantifyingthese processes correspond to the Young's modulus E and Poisson's ratioν, the rock toughness K_(1c), the fracturing fluid viscosity μ (assuminga Newtonian fluid), and the leak-off coefficient C_(l), respectively.There is also the issue of the fluid lag λ, the distance between thefront of the fracturing fluid and the crack edge, which brings into theformulation the magnitude of far-field stress σ_(o) (perpendicular tothe fracture plane) and the virgin pore pressure p_(o).

Multiple embodiments of the present invention are described in thisdisclosure. Some embodiments deal with radial hydraulic fractures, andsome other embodiments deal with plane strain (KGD) fractures, and stillother embodiments are general to all types of fractures. Further,different embodiments employ various scalings and various parametricspaces. For purposes of illustration, and not by way of limitation, theremainder of this disclosure is organized by different types ofparametric spaces, and various other organizational breakdowns areprovided within the discussion of the different types of parametricspaces.

I. Embodiments Utilizing a First Parametric Space

A. Radial Fractures

The problem of a radial hydraulic fracture driven by injecting a viscousfluid from a “point”-source, at a constant volumetric rate Q_(o) isschematically shown in FIGS. 1 and 2. Under conditions where the lag isnegligible (λ/R<<1), determining the solution of this problem consistsof finding the aperture w of the fracture, and the net pressure p (thedifference between the fluid pressure p_(f) and the far-field stressσ_(o)) as a function of both the radial coordinate r and time t, as wellas the evolution of the fracture radius R(t). The functions R(t),w(r,t), and p(r,t) depend on the injection rate Q_(o) and on the 4material parameters E′, μ′, K′, and C′ respectively defined as

$\begin{matrix}{{{E^{\prime} = \frac{E}{1 - v^{2}}}\mspace{14mu} u^{\prime} = {12\;\mu}}\mspace{14mu}{K^{\prime} = {4\left( \frac{2}{\pi} \right)^{1/2}K_{Ic}}}{C^{\prime} = {2C_{I}}}} & (1)\end{matrix}$

The three functions R(t), w(r,t), and p(r,t) are determined by solving aset of equations which can be summarized as follows.

Elasticity Equation

$\begin{matrix}{w = {\frac{R}{E^{\prime}}{\int_{0}^{1}{{G\left( {{r/R},s} \right)}{p\left( {{sR},t} \right)}s\ {\mathbb{d}s}}}}} & (2)\end{matrix}$where G is a known elastic kernel. This singular integral equationexpresses the non-local dependence of the fracture width w on the netpressure p.Lubrication Equation

$\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{1}{r}\frac{\partial\;}{\partial r}\left( {{rw}^{3}\frac{\partial p}{\partial r}} \right)}} & (3)\end{matrix}$

This non-linear differential equation governs the flow of viscousincompressible fluid inside the fracture. The function g(r,t) denotesthe rate of fluid leak-off, which evolves according to

$\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(r)}}}} & (4)\end{matrix}$where t_(o)(r) is the exposure time of point r (i.e., the time at whichthe fracture front was at a distance r from the injection point). Theleak-off law (4) is an approximation with the constant C′ lumpingvarious small scale processes (such as displacement of the pore fluid bythe fracturing fluid). In general, (4) can be defended under conditionswhere the leak-off diffusion length is small compared to the fracturelength.Global Volume Balance

$\begin{matrix}{{Q_{o}t} = {{2\pi{\int_{0}^{R}{{wr}\ {\mathbb{d}r}}}} + {2\pi{\int_{0}{r{\int_{0}^{R{(\tau)}}{{g\left( {r,\tau}\  \right)}{\mathbb{d}r}\ {\mathbb{d}\tau}}}}}}}} & (5)\end{matrix}$This equation expresses that the total volume of fluid injected is equalto the sum of the fracture volume and the volume of fluid lost in therock surrounding the fracture.Propagation Criterion

$\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{R - r}}},{{1 - \frac{r}{R}} ⪡ 1}} & (6)\end{matrix}$Within the framework of linear elastic fracture mechanics, this equationembodies the fact that the fracture is always propagating and thatenergy is dissipated continuously in the creation of new surfaces inrock (at a constant rate per unit surface). Note that (6) implies thatw=0 at the tip.Tip Conditions

$\begin{matrix}{{{w^{3}\frac{\partial p}{\partial r}} = 0},\mspace{31mu}{r = R}} & (7)\end{matrix}$This zero fluid flow rate condition (q=0) at the fracture tip isapplicable only if the fluid is completely filling the fracture(including the tip region) or if the lag is negligible at the scale ofthe fracture. Otherwise, the equations have to be altered to account forthe phenomena taking place in the lag zone as discussed below.Furthermore, the lag size λ(t) is unknown, see FIG. 2.

The formulated model for the radial fracture or similar model for aplanar fracture gives a rigorous account for various physical mechanismsgoverning the propagation of hydraulic fractures, however, is based onnumber of assumptions which may not hold for some specific classes offractures. Particularly, the effect of fracturing fluid buoyancy (thedifference between the density of fracturing fluid and the density ofthe host rock) is one of the driving mechanisms of vertical magma dykes(though, inconsequential for the horizontal disk shaped magma fractures)is not considered in this proposal. Other processes which could berelevant for the hydraulic fracture propagation under certain limitedconditions which are not discussed here include a process zone near thefracture tip, fracturing fluid cooling and solidification effects (asrelevant to magma-driven fractures), capillarity effects at the fluidfront in the fracture, and deviations from the one-dimensional leak-offlaw.

1. Propagation Regimes of Finite Fractures

Scaling laws for finite radial fracture driven by fluid injected at aconstant rate are considered next. Similar scaling can be developed forother geometries and boundary conditions. Regimes with negligible fluidlag are differentiated from regimes with non-negligible fluid lag.

a. Regimes with Negligible Fluid Lag.

Propagation of a hydraulic fracture with zero lag is governed by twocompeting dissipative processes associated with fluid viscosity andsolid toughness, respectively, and two competing components of the fluidbalance associated with fluid storage in the fracture and fluid storagein the surrounding rock (leak-off). Consequently, limiting regimes ofpropagation of a fracture can be associated with dominance of one of thetwo dissipative processes and/or dominance of one of the two fluidstorage mechanisms. Thus, four primary asymptotic regimes of hydraulicfracture propagation with zero lag can be identified where one of thetwo dissipative mechanisms and one of the two fluid storage componentsare vanishing: storage-viscosity (M), storage-toughness (K),leak-off-viscosity ({tilde over (M)}), and leak-off-toughness ({tildeover (K)}) dominated regimes. For example, fluid leak-off is negligiblecompared to the fluid storage in the fracture and the energy dissipatedin the flow of viscous fluid in the fracture is negligible compared tothe energy expended in fracturing the rock in thestorage-viscosity-dominated regime (M). The solution in thestorage-viscosity-dominated regime is given by the zero-toughness,zero-leak-off solution (K′=C′=0). As used herein, the letters M (forviscosity) and K (for toughness) are used to identify which dissipativeprocess is dominant and the symbol tilde (˜) (for leak-off) and no-tilde(for storage in the fracture) are used to identify which fluid balancemechanism is dominant.

Consider general scaling of the finite fracture which hinges on definingthe dimensionless crack opening Ω, net pressure Π, and fracture radius γas:w=εLΩ(ρ;P ₁ ,P ₂), p=εE′Π(ρ;P ₁ ,P ₂), R=γ(P ₁ ,P ₂)L  (8)These definitions introduce a scaled coordinate ρ=r/R(t) (0≦ρ≦1), asmall number ε(t), a length scale L(t) of the same order of magnitude asthe fracture length R(t), and two dimensionless evolution parametersP₁(t) and P₂(t), which depend monotonically on t. The form of thescaling (8) can be motivated from elementary elasticity considerations,by noting that the average aperture scaled by the fracture radius is ofthe same order as the average net pressure scaled by the elasticmodulus.

Four different scalings can be defined to emphasize above differentprimary limiting cases. These scalings yield power law dependence of L,ε, P₁, and P₂ on time t; i.e. L˜t^(α), ε˜t^(δ), P₁˜t^(β) ¹ , P₂˜t^(β) ², see Table 1 for the case of a radial fracture. Furthermore, theevolution parameters can take either the meaning of a toughness (K_(m),K_({tilde over (m)})), or a viscosity (M_(k), M _(k) ), or a storage(S_({tilde over (m)}), S _(k) ), or a leak-off coefficient (C_(m),C_(k)).

TABLE 1 Small parameter ε, lengthscale L, and parameters P₁ and P₂ forthe two storage scalings (viscosity and toughness) and the two leak-offscalings (viscosity and toughness). Scaling ε L P₁ P₂ storage/viscosity(M) $\left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}$$\left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/9}$$K_{m} = {K^{\prime}\left( \frac{t^{2}}{\mu^{\prime 5}E^{\prime 13}Q_{o}^{3}} \right)}^{1/18}$$C_{m} = {\text{C}^{\prime}\left( \frac{E^{\prime 4}t^{7}}{\mu^{\prime 4}Q_{o}^{6}} \right)^{1/18}}$storage/toughness (K)$\left( \frac{K^{\prime 6}}{E^{\prime 6}Q_{o}t} \right)^{1/5}$$\left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/5}$$M_{k} = {\mu^{\prime}\left( \frac{Q_{o}^{3}E^{\prime 13}}{K^{\prime 18}t^{2}} \right)}^{1/5}$$C_{k} = {\text{C}^{\prime}\left( \frac{E^{\prime 8}t^{3}}{K^{\prime 8}Q_{o}^{2}} \right)^{1/10}}$leak-off/viscosity (M)$\left( \frac{\mu^{\prime 4}C^{\prime 6}}{E^{\prime 4}Q_{o}^{2}t^{3}} \right)^{1/16}$$\left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/4}$$K_{\overset{\sim}{m}} = {K^{\prime}\left( \frac{t}{E^{\prime 12}\mu^{\prime 4}C^{\prime 2}Q_{o}^{2}} \right)}^{\frac{1}{16}}$$S_{\overset{\sim}{m}} = \left( \frac{\mu^{\prime 4}Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}t^{7}} \right)^{1/16}$leak-off/toughness (K)$\left( \frac{K^{\prime 8}C^{\prime 2}}{E^{\prime 8}Q_{o}^{2}t} \right)^{1/8}$$\left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/4}$$M_{\overset{\sim}{k}} = {\mu^{\prime}\left( \frac{E^{\prime 12}C^{\prime 2}Q_{o}^{2}}{K^{\prime 16}t} \right)}^{1/4}$$S_{\overset{\sim}{k}} = \left( \frac{K^{\prime 8}Q_{o}^{2}}{E^{\prime 8}C^{\prime 10}t^{3}} \right)^{1/8}$

The regimes of solutions can be conceptualized in a rectangularparametric space MK{tilde over (K)}{tilde over (M)} shown in FIG. 3.Each of the four primary regimes (M, K, {tilde over (M)}, and {tildeover (K)}) of hydraulic fracture propagation corresponding to thevertices of the diagram is dominated by only one component of fluidglobal balance while the other can be neglected (i.e. respective P₁=0,see Table 1) and only one dissipative process while the other can beneglected (i.e. respective P₂=0, see Table 1). The solution for each ofthe primary regimes has the property that it evolves with time taccording to a power law. In particular, the fracture radius R evolvesin these regimes according to l˜t^(α) where the exponent α depends onthe regime of propagation: α= 4/9, ⅖, ¼, ¼ in the M-, K-, {tilde over(M)}-, {tilde over (K)}-regime, respectively. As follows from thestationary tip solution (see below), the behavior of the solution at thetip also depends on the regime of solution: Ω˜(1−ρ)^(2/3) at theM-vertex, Ω˜(1−ρ)^(5/8) at the {tilde over (M)}-vertex, andΩ˜(1−ρ)^(1/2) at the K- and {tilde over (K)}-vertices.

The edges of the rectangular phase diagram MK{tilde over (K)}{tilde over(M)} can be identified with the four secondary limiting regimescorresponding to either the dominance of one of the two fluid globalbalance mechanisms or the dominance of one of the two energy dissipationmechanisms: storage-edge (MK, C_(m)=C_(k)=0), leak-off-edge ({tilde over(M)}{tilde over (K)}, S_({tilde over (m)})=S_({tilde over (k)})=0),viscosity-edge (M{tilde over (M)}, K_(m)=K_({tilde over (m)})=0), andK{tilde over (K)}-toughness-edge (M_(k)=M_(k)=0).

The regime of propagation evolves with time, since the parameters M's,K's, C's and S's depend on t. With respect to the evolution of thesolution in time, it is useful to locate the position of the state pointin the MK{tilde over (K)}{tilde over (M)} space in terms of thedimensionless times τ_(mk)=t/t_(mk),τ_({tilde over (m)}{tilde over (k)})=t/t_({tilde over (m)}{tilde over (k)}),τ_({tilde over (m)}{tilde over (m)})=t/t_({tilde over (m)}{tilde over (m)}),and τ_(k k) =t/t_(k k) where the time scales are defined as

$\begin{matrix}{{t_{mk} = \left( \frac{\mu^{\prime 5}E^{\prime 13}Q_{o}^{3}}{K^{\prime 18}} \right)^{1/2}},} & \left( {9a} \right) \\{t_{\overset{\sim}{m}\overset{\sim}{k}} = \frac{\mu^{\prime 4}E^{\prime 2}Q_{o}^{2}}{K^{\prime 6}}} & \left( {9b} \right) \\{t_{m\overset{\sim}{m}} = \left( \frac{\mu^{\prime 4}Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}} \right)^{1/7}} & \left( {9c} \right) \\{t_{k\overset{\sim}{k}} = \left( \frac{K^{\prime 8}Q_{o}^{2}}{E^{\prime 8}C^{\prime 10}} \right)^{1/3}} & \left( {9d} \right)\end{matrix}$Only two of these times are independent, however, sincet_({tilde over (m)}{tilde over (k)})=t_(mk) ^(8/5)t_(k{tilde over (k)})^(−3/5) and t_(m{tilde over (m)})=t_(mk) ^(8/35)t_(k{tilde over (k)})^(27/35). Note that the parameters M's, K's, C's and S's can be simplyexpressed in terms of these times according toK_(m)=M_(k) ^(−5/18)=τ_(mk) ^(1/9),K_({tilde over (m)})=M_({tilde over (k)})^(−1/4)=τ_({tilde over (m)}{tilde over (k)}) ^(1/16),C_(m)=S_({tilde over (m)}) ^(−8/9)=τ_(m{tilde over (m)}) ^(7/18),C_(k)=S_({tilde over (k)}) ^(−4/5)=τ_(k{tilde over (k)}) ^(3/10)  (10)

The dimensionless times τ'S define evolution of the solution along therespective edges of the rectangular space MK{tilde over (K)}{tilde over(M)}. A point in the parametric space MK{tilde over (K)}{tilde over (M)}is thus completely defined by any pair combination of these four times,say (τ_(mk), τ_(k k) ). The position (τ_(mk), τ_(k k) ) of the statepoint can in fact be conceptualized at the intersection of two rays,perpendicular to the storage- and toughness-edges respectively.Furthermore, the evolution of the solution regime in the MK{tilde over(K)}{tilde over (M)} space takes place along a trajectory correspondingto a constant value of the parameter η, which is related to the ratiosof characteristic times

$\begin{matrix}{{\eta = \frac{E^{{\prime 11}/2}\mu^{{\prime 3}/2}C^{\prime 2}Q_{o}^{1/2}}{K^{\prime 7}}},{{{with}\mspace{14mu}\frac{t_{m\overset{\sim}{m}}}{t_{mk}}} = \eta},{\frac{t_{m\overset{\sim}{m}}}{t_{mk}} = \eta},{\frac{t_{k\overset{\sim}{k}}}{t_{mk}} = \eta^{{- 5}/3}},{\frac{t_{m\overset{\sim}{m}}}{t_{k\overset{\sim}{k}}} = \eta^{8/21}}} & (11)\end{matrix}$(One of Such Trajectories is Shown at 310 in FIG. 3).

In view of the dependence of the parameters M's, K's, C's, and S's ontime, (10), the M-vertex corresponds to the origin of time, and the{tilde over (K)}-vertex to the end of time (except for an impermeablerock). Thus, given all the problem parameters which completely definethe number η, the system evolves with time (say time τ_(mk)) along aη-trajectory, starting from the M-vertex (K_(m)=0, C_(m)=0) and endingat the {tilde over (K)}-vertex (M _(k) =0, S _(k) =0). If η=0, twopossibilities exist: either the rock is impermeable (C′=0) and thesystem evolves along the storage edge from M to K, or the fluid isinviscid (μ′=0) and the system then evolves along the toughness-edgefrom K to {tilde over (K)}. If η=∞, then either K′=0 (corresponding to apre-existing discontinuity), and the system evolves along theviscosity-edge from M to {tilde over (M)}; or C′=∞ (corresponding tozero fluid storage in the fracture) and the system evolves along theleak-off-edge from the {tilde over (M)} to the {tilde over (K)}. Thuswhen η is decreasing (which can be interpreted for example as andecreasing ratio t_(m{tilde over (m)})/t_(mk)), the trajectory isattracted by the K-vertex, and when η is increasing the trajectory isattracted to the {tilde over (M)}-vertex. The dependence of the scaledsolution F can thus be expressed in the form F(ρ,τ;η), where τ is one ofthe dimensionless time, irrespective of the adopted scaling.

b. Regimes with Non-negligible Fluid Lag.

Under certain conditions (e.g., when a fracture propagates alongpre-existing discontinuity K′=0 and confining stress σ_(o), is smallenough), the length of the lag between the crack tip and the fluid frontcannot be neglected with respect to the fracture size. In someembodiments of the present invention, fluid pressure in the lag zone canbe considered to be zero compared to the far-field stress σ_(o), eitherbecause the rock is impermeable or because there is cavitation of thepore fluid. Under these conditions, the presence of the lag brings σ_(o)in the problem description, through an additional evolution parameterP₃(t), which is denoted T_(m) in the M-scaling (or T_({tilde over (m)})in the {tilde over (M)}-scaling) and has the meaning of dimensionlessconfining stress. This extra parameter can be expressed in terms of anadditional dimensionless time as

$\begin{matrix}{{T_{m} = \tau_{om}^{1/3}}{{{with}\mspace{14mu}\tau_{om}} = {{\frac{t}{t_{om}}\mspace{14mu}{and}\mspace{14mu} t_{om}} = \frac{\mu^{\prime}E^{\prime 2}}{\sigma_{o}^{3}}}}} & (12)\end{matrix}$Now the parametric space can be envisioned as the pyramid MK{tilde over(K)}{tilde over (M)}-OÕ, depicted in FIG. 4, with the position of thestate point identified by a triplet, e.g., (T_(m), K_(m), C_(k)) or(τ_(om), τ_(mk), τ_(k k) ). In accord with the discussion of the zerolag case, OÕ-edge corresponds to the viscosity-dominated regime(K_(m)=K_({tilde over (m)})=0) under condition of vanishing confiningstress (T_(m)=T_({tilde over (m)})=0), where the endpoints, O- andÕ-vertices correspond to the limits of storage and leak-off-dominatedcases.

The system evolves from the O-vertex towards the {tilde over (K)}-vertexfollowing a trajectory which depends on all the parameters of theproblem (410, FIG. 4). The trajectory depends on two numbers which canbe taken as η defined in (11) (independent of σ_(o)) andφ=t_(om)/t_(mk). It should be noted that the O-vertex from wherefracture evolution initiates is a singular point as (i) it correspondsto the infinitely fast initial fracture propagation (propagation of anunconfined fracture, σ_(o)=0, along preexisting discontinuity, K′=0)(ii) it corresponds to the infinite multitude of self-similar solutionsparameterized by the ray along which the solution trajectory is emergingfrom the O-vertex.

If φ<<1 and φ<<η (e.g. the confining stress σ_(o) is “large”), thetrajectory follows essentially the OM-edge, and then from the M-vertexremains within the MK{tilde over (K)}{tilde over (M)}-rectangle.Furthermore, the transition from O to M takes place extremely morerapidly than the evolution from the M to the {tilde over (K)}-vertexalong a η-trajectory (or from M to the K-vertex if the rock isimpermeable). In other words, the parametric space can be reduced to theMK{tilde over (K)}{tilde over (M)}-rectangle, and the lag can thus beneglected if φ<<1 and φ<<η. Through this reduction in the dimensions ofthe parametric space, the M-vertex becomes the apparent starting pointof the evolution of a fluid-driven fracture without lag. The “penalty”for this reduction is a multiple boundary layer structure of thesolution near the M-vertex.

If the rock is impermeable (C′=0), the solution is restricted to evolveon the MKO face of the parametric space (see FIG. 5), from O to Kfollowing a φ-trajectory 510. However, there is no additional time scaleassociated with the OK-edge and thus the transition OK takes place“rapidly” if φ>>1; this is a limiting case where the lag can beneglected, as the solution is always in the asymptotic K-regime.

2. Structure of the Solution Near the Tip of Propagating HydraulicFracture

The nature of the solution near the tip of a propagating fluid-drivenfracture can be investigated by analyzing the problem of a semi-infinitefracture propagating at a constant speed V, see FIGS. 6 and 7. In thefollowing, a distinction is made between regimes/scales with negligibleand non-negligible lag between the crack tip and the fluid front.Although a lag of a prior unknown length λ between the crack tip and thefluid front must necessarily exist on a physical ground, as otherwisethe fluid pressure at the tip has negative singularity, there arecircumstances where λ is small enough compared to the relevantlengthscales that it can be neglected. (This issue is similar to the useof the solutions of linear elastic fracture mechanics which yield“unphysical” stress singularity at the fracture tip). In theseregimes/scales, the solution is characterized by a singular behavior,with the nature of the singularity being a function of the problemparameters and the scale of reference.

a. Regimes/Scales with Negligible Fluid Lag.

In view of the stationary nature of the considered tip problem, thefracture opening ŵ, net pressure {circumflex over (p)} and flow rate{circumflex over (q)} are only a function of the moving coordinate{circumflex over (x)}, see FIGS. 6 and 7. The system of equationsgoverning ŵ({circumflex over (x)}) and {circumflex over (p)}({circumflexover (x)}) can be written as

$\begin{matrix}{{\hat{p} = {\frac{E^{\prime}}{4\;\pi}{\int_{0}^{\infty}{\frac{\mathbb{d}\hat{w}}{\mathbb{d}\hat{s}}\frac{d\hat{s}}{\hat{x} - \hat{s}}}}}}\ ,{{\frac{{\hat{w}}^{3}}{\mu^{\prime}}\frac{\mathbb{d}\hat{p}}{\mathbb{d}\hat{x}}} = {{V\hat{w}} + {2C^{\prime}V^{1/2}{\hat{x}}^{1/2}}}},{{\lim\limits_{\hat{x}\rightarrow 0}\frac{\hat{w}}{{\hat{x}}^{1/2}}} = \frac{K^{\prime}}{E^{\prime}}},{{\lim\limits_{\hat{x}\rightarrow 0}{{\hat{w}}^{3}\frac{\mathbb{d}\hat{p}}{\mathbb{d}\hat{x}}}} = 0}} & (13)\end{matrix}$

The singular integral equation (13)_(a) derives from elasticity, whilethe Reynolds equation (13)_(b) is deduced from the Poiseuille({circumflex over (q)}=ŵ³/μ′d{circumflex over (p)}/d{circumflex over(x)}), continuity (V dŵ/d{circumflex over (x)}−d{circumflex over(q)}/d{circumflex over (x)}+ĝ=0), and Carter's leak-off laws(ĝ=C′√{square root over (V/{circumflex over (x)})}). Equation (13)_(c)expresses the crack propagation criterion, while the zero flow ratecondition at the tip, (13)_(d), arises from the assumption of zero lag.

Analogously to the considerations for the finite fracture, four primarylimiting regimes of propagation of a semi-infinite fracture with zerolag can be identified where one of the two dissipative mechanisms andone of the two fluid storage components are vanishing: storage-viscosity(m), storage-toughness (k), leak-off-viscosity ({tilde over (m)}), andleak-off-toughness ({tilde over (k)}) dominated regimes. Each of theregimes correspond to the respective vertex of the rectangularparametric space of the semi-infinite fracture. However, in the contextof the semi-infinite fracture, the storage-toughness (k) andleak-off-toughness ({tilde over (k)}) dominated regimes are identicalsince the corresponding zero viscosity (μ′=0) solution of (13) isindependent of the balance between the fluid storage and leak-off, andis given by the classical linear elastic fracture mechanics (LEFM)solution ŵ=(K′/E′){circumflex over (x)}^(1/2) and {circumflex over(p)}=0. Therefore, the toughness edge k{tilde over (k)} of therectangular parameteric space for the semi-infinite fracture collapsesinto a point, which can be identified with either k- or {tilde over(k)}-vertex, and the rectangular space itself into the triangularparametric space mk{tilde over (m)}, see FIG. 7.

The primary storage-viscosity, toughness, and leak-off-viscosityscalings associated with the three primary limiting regimes (m, k or{tilde over (k)}, and {tilde over (m)}) are as follows

$\begin{matrix}{{{\hat{\xi}}_{m,k,\overset{\sim}{m}} = \frac{\hat{x}}{l_{m,k,\overset{\sim}{m}}}},{{\hat{\Omega}}_{m,k,\overset{\sim}{m}} = \frac{\hat{w}}{l_{m,k,\overset{\sim}{m}}}},{{\hat{\Pi}}_{m,k,\overset{\sim}{m}}\; = \frac{\hat{p}}{E^{\prime}}},{{\hat{\Psi}}_{m,k,\overset{\sim}{m}} = \frac{\hat{q}}{V\; l_{m,k,\overset{\sim}{m}}}}} & (14)\end{matrix}$where the three lengthscales l_(m), l_(k), and l_({tilde over (m)}), aredefined as l_(m)=μ′V/E′, l_(k)=(K′/E′)²,l_({tilde over (m)})=V^(1/3)(2μ′C′)^(2/3)/E′^(2/3). The solutionF={{circumflex over (Ω)},{circumflex over (Π)}} in the various scalingscan be shown to be of the form {circumflex over (F)}_(m)({circumflexover (ξ)}_(m);c_(m),k_(m)), {circumflex over (F)}_(k)({circumflex over(ξ)}_(k);m_(k),m _(k) ), {circumflex over(F)}_({tilde over (m)})(ξ_({tilde over (m)});s_({tilde over (m)}),k_({tilde over (m)})),with the letters m's, k's, s's and c's representing dimensionlessviscosity, toughness, storage, and leak-off coefficient, respectively.k _(m) =m _(k) ^(−1/2)=(l _(k) l _(m))^(1/2) , k _({tilde over (m)}) =m_({tilde over (k)}) ^(−1/3)=(l _(k) l _({tilde over (m)}))^(1/2) , c_(m) =s _({tilde over (m)}) ^(−3/2)=(l _({tilde over (m)}) l_(m))^(3/2)  (15)

For example, a point in the mk{tilde over (m)} ternary diagramcorresponds to a certain pair (k_(m), c_(m)) in the viscosity scaling,with the m-vertex corresponding to c_(m)=0 and k_(m)=0. The vertexsolutions (denoted by the subscript ‘0’) are given by{circumflex over (Ω)}_(m0)=β_(m0){circumflex over (ξ)}_(m) ^(2/3),{circumflex over (Π)}_(m0)=δ_(m0){circumflex over (ξ)}_(m) ^(−1/3);{circumflex over (Ω)}_(k0)={circumflex over (ξ)}_(k) ^(1/2), {circumflexover (Π)}_(k0)=0; and{circumflex over(Ω)}_({tilde over (m)}0)=β_({tilde over (m)}0){circumflex over(ξ)}_({tilde over (m)}) ^(5/8), {circumflex over(Π)}_({tilde over (m)}0)=δ_({tilde over (m)}0){circumflex over(ξ)}_({tilde over (m)}) ^(−3/8)  (16)with β_(m0)=2^(1/3)3^(5/6), δ_(m0)=−6^(−2/3),β_({tilde over (m)}0)≅2.534, δ_({tilde over (m)}0)≅−0.164. Thus whenthere is only viscous dissipation (edge m{tilde over (m)} correspondingto fracture propagation along preexisting discontiuity K′=0) the tipbehavior is of the form ŵ˜{circumflex over (x)}^(2/3), {circumflex over(p)}˜−{circumflex over (x)}^(−1/3) in the storage-dominated case,m-vertex, (impermeable rock C′=0) and of the form ŵ˜{circumflex over(x)}^(5/8), {circumflex over (p)}˜−{circumflex over (x)}^(−3/8) in theleak-off dominated case, {tilde over (m)}-vertex. On the other hand, thek-vertex pertains to a fracture driven by an inviscid fluid (μ′=0); thisvertex is associated with the classical tip solution of linear elasticfracture mechanics ŵ˜{circumflex over (x)}^(1/2). The general case of afluid-driven fracture with no leak-off (C′=0) or negligible storagenaturally corresponds to the mk- or {tilde over (m)}k-edges,respectively. However, a more general interpretation of the mk{tildeover (m)} parametric space can be seen by expressing the numbers m's,k's, s's, and c's in terms of a dimensionless velocity ν, and aparameter η which only depends on the parameters characterizing thesolid and the fluid

$\begin{matrix}{{v = {\frac{l_{m}}{l_{k}} = \frac{V}{V*}}},{\hat{\eta} = {\frac{l_{m}l_{k}^{2}}{l_{\overset{\sim}{m}}^{3}} = \frac{4E^{\prime 3}\mu^{\prime}C^{\prime 2}}{K^{\prime 4}}}}} & (17)\end{matrix}$where V*=K′²/μ′E is a characteristic velocity. Hence, k_(m)=ν^(−1/2),k_({tilde over (m)})={circumflex over (η)}^(−1/6)ν^(−1/6),c_(m)={circumflex over (η)}^(1/2)ν⁻¹. The above expressions indicatethat the solution moves from the m-vertex towards the k-vertex withdecreasing dimensionless velocity ν, along a trajectory which dependsonly {circumflex over (η)}. With increasing {circumflex over (η)}, thetrajectory is pulled towards the {tilde over (m)}-vertex. Since the tipvelocity of a finite fracture decreases with time (at least underconstant injection rate), the tip solution interpreted from thisstationary solution is seen to evolve with time. In other words, as thelength scales l_(m) and l_({tilde over (m)}) evolve with time, thenature of the solution in the tip region at a given physical scaleevolves accordingly.

The solution along the edges of the mkm-triangle, namely, the viscositymm-edge (k_(m)=0 or k_({tilde over (m)})=0), the storage mk-edge(c_(m)=0 or m _(k) =0), and the {tilde over (m)}k-edge(s_({tilde over (m)})=0 or m_(k)=0) has been obtained both in the formof series expansion in the neighborhood of the vertices and numericallyfor finite values of the non-zero parameters. These results wereobtained in part by recognizing that the solution can be furtherrescaled along each edge to eliminate the remaining parameter. Forexample, the tip solution along the mk-edge, which is governed byparameter k_(m) in the m-scaling, upon rescaling to the mixed scalingcan be expressed as {circumflex over (F)}_(mk)({circumflex over(ξ)}_(mk)) where {circumflex over (ξ)}_(mk)={circumflex over (x)}/l_(mk)with l_(mk)=l_(k) ³/l_(m) ².

The m{tilde over (m)}-, mk-, and {tilde over (m)}k-solutions obtained sofar give a glimpse on the changing structure of the tip solution atvarious scales, and how these scales change with the problem parameters,in particular with the tip velocity ν. Consider for example themk-solution (edge of the triangle corresponding to the case ofimpermeable rock) for the opening {circumflex over (Ω)}_(mk)({circumflexover (ξ)}_(mk)), with {circumflex over (Ω)}_(mk)=k_(m) ⁻⁴{circumflexover (Ω)}_(m)=m_(k){circumflex over (Ω)}_(k). Expansion of the{circumflex over (Ω)}_(mk) at {circumflex over (ξ)}_(mk)=0 and at{circumflex over (ξ)}_(mk)=∞ is of the form{circumflex over (Ω)}_(mk)=β_(m0){circumflex over (ξ)}_(mk)^(2/3)+β_(mk1){circumflex over (ξ)}_(mk) ^(h) +O({circumflex over(ξ)}_(mk) ^(−1/3)) at{circumflex over (ξ)}_(mk)=∞, {circumflex over(Ω)}_(mk={circumflex over (ξ)}) _(mk) ^(1/2)+β_(km1){circumflex over(ξ)}_(mk) +O({circumflex over (ξ)}_(mk) ^(3/2))  (18)

The exponent h≅0.139 in the “alien” term {circumflex over (ξ)}_(mk) ^(h)of the far-field expansion (18)₁ is the solution of certaintranscendental equation obtained in connection with correspondingboundary layer structure. In this case, the boundary layer arisesbecause ŵ˜{circumflex over (x)}^(1/2) near {circumflex over (x)}=0 ifK′>0, but ŵ˜{circumflex over (x)}^(2/3) when K′=0. The behavior of themk-solution at infinity corresponds to the m-vertex solution. Themk-solution shows that {circumflex over (Ω)}_(mk)≅β_(m0){circumflex over(ξ)}_(m) ^(2/3) for {circumflex over (ξ)}_(mk)>{circumflex over(ξ)}_(mk∞), with {circumflex over (ξ)}_(mk∞)=O(1), with the consequencethat there will be corresponding practical range of parameters for whichthe global solution for C′=0 is characterized by the m-vertex asymptoticbehavior ŵ˜{circumflex over (x)}^(2/3), {circumflex over(p)}˜−{circumflex over (x)}^(−1/3) (viscous dissipation only), althoughŵ˜{circumflex over (x)}^(1/2) in a very small region near the tip.Taking for example V≅1 m/s, E≅10³ MPa, μ′≅10⁻⁶ MPa·s, K′≅1 MPa·m^(1/2),and C′=0, then l_(mk)≅10⁻² m. Hence, at distance larger than 10⁻² m, thesolution behaves as if the impermeable rock has no toughness and thereis only viscous dissipation. As discussed further below, the m-vertexsolution develops as an intermediate asymptote at some small distancefrom the tip in the finite fracture, provided the lengthscale l_(mk) ismuch smaller than the fracture dimension R.

b. Regimes/Scales with Non-negligible Fluid Lag.

The stationary problem of a semi-infinite crack propagating at constantvelocity V is now considered, taking into consideration the existence ofa lag of a priori unknown length λ between the crack tip and the fluidfront, see FIG. 2. First, considerations are restricted to impermeablerocks. In this case, the tip cavity is filled with fluid vapors, whichcan be assumed to be at zero pressure.

This problem benefits from different scalings in part because thefar-field stress σ_(o) directly influences the solution, through thelag. Consider for example the mixed stress/storage/viscosity scaling(om)

$\begin{matrix}{{{\hat{\xi}}_{om} = {\frac{\hat{x}}{l_{om}} = {ɛ^{3}{\hat{\xi}}_{m}}}},{{\hat{\Omega}}_{om} = {ɛ^{2}{\hat{\Omega}}_{m}}},{{\hat{\Pi}}_{om} = {ɛ^{- 1}{\hat{\Pi}}_{m}}},{{{with}\mspace{14mu} l_{om}} = {ɛ^{- 3}l_{m}}},{ɛ = \frac{\sigma_{o}}{E^{\prime}}}} & (19)\end{matrix}$

It can be shown that the solution is of the form {circumflex over(F)}_(om)({circumflex over (ξ)}_(om);k_(om)), where k_(om)=ε^(1/2) k_(m)is the dimensionless toughness in this new scaling; {circumflex over(F)}_(om) behaves according to the k-vertex asymptote near the tip({circumflex over (Ω)}_(om)≅k_(om){circumflex over (ξ)}_(om) ^(1/2) near{circumflex over (ξ)}_(om)=0) and to the m-vertex asymptote far awayfrom the tip ({circumflex over (Ω)}_(om)≅β_(om){circumflex over(ξ)}_(om) ^(2/3) for {circumflex over (ξ)}_(om)>>1). The scaled lagλ_(om)=λ/l_(om) continuously decreases with k_(om), from a maximum valueλ_(mso)=0.36 reached either when K′=0 or σ_(o)=0. The decrease of λ_(om)with k_(om) becomes exponentially fast for large toughness (practicallywhen k_(om)≳4). Furthermore, analysis of the solution indicates that{circumflex over (F)}_(om)({circumflex over (ξ)}_(om);k_(om)) can berescaled into {circumflex over (F)}_(mk)({circumflex over (ξ)}_(mk)) forlarge toughness (k_(om)≳4){circumflex over (ξ)}_(mk)=k_(om) ⁻⁶{circumflex over (ξ)}_(om),{circumflex over (Ω)}_(mk)=k_(om) ⁻⁴{circumflex over (Ω)}_(om),{circumflex over (Π)}_(mk)=k_(om) ²{circumflex over (Π)}_(mk)  (20)

These considerations show that within the context of the stationary tipsolution the fluid lag becomes irrelevant at the scales of interest ifk_(om)≳4, and can thus be assumed to be zero (with the implication that{circumflex over (q)}=0 at the tip, which leads to a singularity of thefluid pressure.) Also, the solution becomes independent of the far-fieldstress σ_(o) when k_(om)≳4 (except as a reference value of the fluidpressure) and it can be mapped within the mk{tilde over (m)} parametricspace introduced earlier.

In permeable rocks, pore fluid is exchanged between the tip cavity andthe porous rock and flow of pore fluid within the cavity is takingplace. The fluid pressure in the tip cavity is thus unknown andfurthermore not uniform. Indeed, pore fluid is drawn in by suction atthe tip of the advancing fracture, and is reinjected to the porousmedium behind the tip, near the interface between the two fluids. (Porefluid must necessarily be returning to the porous rock from the cavity,as it would otherwise cause an increase of the lag between thefracturing fluid and the tip of the fracture, and would thus eventuallycause the fracture to stop propagating). Only elements of the solutionfor this problem exists so far, in the form of a detailed analysis ofthe tip cavity under the assumption that ŵ˜{circumflex over (x)}^(1/2)in the cavity.

This analysis shows that the fluid pressure in the lag zone can beexpressed in terms of two parameters: a dimensionless fracture velocityν=Vλ/c and a dimensionless rock permeability ç=kE′³/(λ^(1/2)K′³), wherek and c denote respectively the intrinsic rock permeability anddiffusivity. Furthermore, the solution is bounded by two asymptoticregimes: drained with the fluid pressure in the lag equilibrated withthe ambient pore pressure p_(o) ( ν<<1 and ç>>1), and undrained with thefluid pressure corresponding to its instantaneous (undrained) value atthe moving fracture tip

$\begin{matrix}{p_{f{({tip})}} = {p_{o} - {\frac{1}{2}\frac{K^{\prime}}{E^{\prime}}\frac{\mu_{o}}{k}\sqrt{\pi\; c\; V}}}} & (21)\end{matrix}$where μ_(o) is the viscosity of the pore fluid. The above expression forP_(f(tip)) indicates that pore fluid cavitation can take place in thelag. Analysis of the regimes of solution suggests that the pore fluidpressure in the lag zone drop below cavitation limit in a wide range ofparameters relevant for propagation of hydraulic fractures and magmadykes, implying a net-pressure lag condition identical to the one forimpermeable rock. This condition allows one to envision the parametricspace for the tip problem in the general case of the permeable rock(leak-off) and the lag (finiteness of the confining stress) as thepyramid mk{tilde over (m)}-oõ, where similarly to the case of the finitefracture, see FIG. 4, vertices o- and õ-correspond to the limits ofstorage and leak-off dominated cases under conditions of vanishingtoughness and confining stress. The stationary tip solution near the om-and õ{tilde over (m)}-edges behaves as k-vertex asymptote (ŵ˜{circumflexover (x)}^(1/2)) near the tip and as the m-vertex (ŵ˜{circumflex over(x)}^(2/3)) and m-vertex (ŵ˜{circumflex over (x)}^(5/8)) asymptote,respectively, far away from the tip.3. Local Tip and Global Structure of the Solution

The development of the general solution corresponding to the arbitraryη-trajectory in the MK{tilde over (K)}{tilde over (M)} rectangle (or(η,φ)-trajectory in the MK{tilde over (K)}{tilde over (M)}-OÕ pyramid)is aided by understanding the asymptotic behavior of the solution in thevicinities of the rectangle (pyramid) vertices and edges. Theseasymptotic solutions can be obtained (semi-) analytically via regular orsingular perturbation analysis. Construction of those solutions to thenext order in the small parameter(s) associated with the respective edge(or vertex) can identify the physically meaningful range of parametersfor which the fluid-driven fracture propagates in the respectiveasymptotic regime (and thus can be approximated by the respective edge(vertex) asymptotic solution). Since the solution trajectory evolveswith time from M-vertex to the {tilde over (K)}-vertex inside of theMK{tilde over (K)}{tilde over (M)}-rectangle (or generally, from theO-vertex to the {tilde over (K)}-vertex inside of the MK{tilde over(K)}{tilde over (M)}-OÕ pyramid), it is helpful to have valid asymptoticsolutions developed in the vicinities of these vertices. The solution inthe vicinity of the some of the vertices (e.g., O, K, and {tilde over(K)}) is a regular perturbation problem, which has been solved for theK-vertex along the MK- and KO-edge of the pyramid. The solution in thevicinity of the M-vertex is challenging since it constitutes a singularperturbation problem for a system of non-linear, non-local equations inmore than one small parameter, namely, K_(m) (along the MK-edge), C_(m)(along the MM-edge), and, generally, E_(m)=T_(m) ⁻¹ (along the MO-edge),given that the nature of the tip singularity changes with a smallperturbation from zero in any of these parameters. Indeed, solution atM-vertex is given by the zero-toughness (K_(m)=0), zero leak-off(C_(m)=0), zero-lag (E_(m)=T_(m) ⁻¹=0) solution which near tip behavioris given by the m-vertex tip solution, Ω_(m)˜(1−ρ)^(2/3) andΠ_(m)˜−(1−ρ)^(−1/3). Small perturbation of the M-vertex in eithertoughness K_(m), or leak-off C_(m), or lag E_(m) changes the nature ofthe near tip behavior to either the toughness asymptoteΩ_(m)˜(1−ρ)^(1/2), or the leak-off asymptote Ω_(m)˜(1−ρ)^(5/8), or thelag asymptote Ω_(m)˜(1−ρ)^(3/2), respectively. This indicates theemergence of the near tip boundary layer (BL) which incorporates arisingtoughness singularity and/or leak-off singularity and/or the fluid lag.If the perturbation is small enough, there exists a lengthscaleintermediate to the fracture length and the BL “thickness” where theouter solution (i.e. the solution away from the fracture tip) and the BLsolution (given by the stationary tip solution discussed above) can bematched to form the composite solution uniformly valid along thefracture. Matching requires that the asymptotic expansions of the outerand the BL solutions over the intermediate lengthscale are identical.

As an illustration, the non-trivial structure of the global solution inthe vicinity of the M-vertex along the MK-edge (i.e., singularperturbation problem in K_(m), while C_(m)=E_(m)=0) is now outlined,corresponding to the case of a fracture in impermeable rock and largeconfining stress (or time). The outer expansion for Ω, Π, anddimensionless fracture radius γ are perturbation expansions in powers ofK_(m) ^(b), b>0. Here the matching not only gives the coefficients inthe expansion, but also determines the exponent b. It can be shown thatthe tip solution along the mk-edge (18) corresponds to the O(1) term inthe inner (boundary layer) expansion at the tip. The inner and outer(global) scaling for the radial fracture are related as

$\begin{matrix}{{{1 - \rho} = {\frac{81}{16\;\gamma_{m\; 0}^{3}}K_{m}^{6}{\hat{\xi}}_{mk}}},{\Omega_{m} = {\frac{9}{4\;\gamma_{m\; 0}}K_{m}^{4}{\hat{\Omega}}_{mk}}},{\Pi_{m} = {\frac{4\;\gamma_{m\; 0}}{9}K_{m}^{- 2}{\hat{\Pi}}_{mk}}}} & (22)\end{matrix}$where γ_(m0) is the O(1) term of the outer expansion for γ given by theM-vertex solution (K_(m)=C_(m)=E_(m)=0). Using the asymptotic expression(18) together with the scaling (22), one finds that the outer and innersolutions match under the condition K_(m) ⁶<<1. Then the leading orderinner and outer solutions form a single composite solution of O(1)uniformly valid along the fracture. That is, to leading order there is alengthscale intermediate to the tip boundary layer thickness K_(m) ⁶Rand the fracture radius R, over which the inner and outer solutionsposses the same intermediate asymptote, corresponding to the m-vertexsolution (16)₁. This solution structure corresponds to the outerzero-toughness solution valid on the lengthscale of the fracture, andthin tip boundary layer given by the mk-edge solution.

To leading order the condition K_(m) ⁶<<1 is merely a condition for theexistence of the boundary layer solution. In order to move away from theM-vertex solution away from the tip, one has to determine the exponent bin the next term in the asymptotic expansion. From this value of b wedetermine the asymptotic validity of the approximation. This can beobtained from the next-order matching between the near tip asymptote inthe outer expansion and the away from tip behavior of the innersolution, see (18). Here the matching to the next order of the outer andinner solutions does not require the next-order inner solution, as thenext order outer solution is matched with the leading order term of theinner solution. The latter appears to be a consequence of the non-localcharacter of the perturbation problem. Then using (18) an expression forthe exponent b=4−6h is obtained which yields b≅3.18 and consequently thenext order contribution in the asymptotic expansion away from the tip.The range of dimensionless toughness in which fracture global (outer)solution can be approximated by the M-vertex solution is, therefore,given by K_(m) ^(3.18)<<1.

B. Plane Strain (KGD) Fractures

The problem of a KGD hydraulic fracture driven by injecting a viscousfluid from a “point”-source, at a constant volumetric rate Q_(o) isschematically shown in FIG. 8. Under conditions where the lag isnegligible, determining the solution of this problem consists of findingthe aperture w of the fracture, and the net pressure p (the differencebetween the fluid pressure p_(f) and the far-field stress σ_(o)) as afunction of both the coordinate x and time t, as well as the evolutionof the fracture radius l(t). The functions l(t), w(x,t), and p(x,t)depend on the injection rate Q_(o) and on the 4 material parameters E′,μ′, K′, and C′ respectively defined as

$\begin{matrix}{E^{\prime} = {{\frac{E}{1 - v^{2}}\mspace{14mu}\mu^{\prime}} = {{12\;\mu\mspace{14mu} K^{\prime}} = {{4\left( \frac{2}{\pi} \right)^{1/2}K_{Ic}\mspace{14mu} C^{\prime}} = {2C_{I}}}}}} & (23)\end{matrix}$The three functions l(t), w(x,t), and p(x,t) are determined by solving aset of equations which can be summarized as follows.Elasticity Equation

$\begin{matrix}{{p\left( {x,t} \right)} = {{- \frac{E^{\prime}}{4\;\pi}}{\int_{- l}^{l}{\frac{\partial{w\left( {s,t} \right)}}{\partial s}\frac{\mathbb{d}s}{s - x}}}}} & (24)\end{matrix}$This singular integral equation expresses the non-local dependence ofthe fracture width w on the net pressure p.Lubrication Equation

$\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{\partial\;}{\partial x}\left( {w^{3}\frac{\partial p}{\partial x}} \right)}} & (25)\end{matrix}$This non-linear differential equation governs the flow of viscousincompressible fluid inside the fracture. The function g(x,t) denotesthe rate of fluid leak-off, which evolves according to

$\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(x)}}}} & (26)\end{matrix}$where t_(o)(x) is the exposure time of point x (i.e., the time at whichthe fracture front was at a distance x from the injection point).Global Volume Balance

$\begin{matrix}{{Q_{o}t} = {{2{\int_{0}^{l}{w\ {\mathbb{d}x}}}} + {2{\int_{0}{\int_{0}^{l{(\tau)}}{{g\left( {x,\tau}\  \right)}{\mathbb{d}x}\ {\mathbb{d}\tau}}}}}}} & (27)\end{matrix}$This equation expresses that the total volume of fluid injected is equalto the sum of the fracture volume and the volume of fluid lost in therock surrounding the fracture.Propagation Criterion

$\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{l - x}}},{{1 - \frac{x}{l}} ⪡ 1}} & (28)\end{matrix}$Within the framework of linear elastic fracture mechanics, this equationembodies the fact that the fracture is always propagating and thatenergy is dissipated continuously in the creation of new surfaces inrock (at a constant rate per unit surface). Note that (28) implies thatw=0 at the tip.Tip Conditions

$\begin{matrix}{{{w^{3}\frac{\partial p}{\partial x}} = 0},{x = l}} & (29)\end{matrix}$This zero fluid flow rate condition (q=0) at the fracture tip isapplicable only if the fluid is completely filling the fracture(including the tip region) or if the lag is negligible at the scale ofthe fracture.1. Propagation Regimes of a KGD Fracture

Propagation of a hydraulic fracture with zero lag is governed by twocompeting dissipative processes associated with fluid viscosity andsolid toughness, respectively, and two competing components of the fluidbalance associated with fluid storage in the fracture and fluid storagein the surrounding rock (leak-off). Consequently, the limiting regimesof propagation of a fracture can be associated with the dominance of oneof the two dissipative processes and/or the dominance of one of the twofluid storage mechanisms. Thus, four primary asymptotic regimes ofhydraulic fracture propagation with zero lag can be identified where oneof the two dissipative mechanisms and one of the two fluid storagecomponents are vanishing: storage-viscosity (M), storage-toughness (K),leak-off-viscosity ({tilde over (M)}), and leak-off-toughness ({tildeover (K)}) dominated regimes. For example, fluid leak-off is negligiblecompared to the fluid storage in the fracture and the energy dissipatedin the flow of viscous fluid in the fracture is negligible compared tothe energy expended in fracturing the rock in thestorage-viscosity-dominated regime (M). The solution in thestorage-viscosity-dominated regime is given by the zero-toughness,zero-leak-off solution (K′=C′=0).

Consider the general scaling of the finite fracture, which hinges ondefining the dimensionless crack opening Ω, net pressure Π, and fractureradius γ asw=εLΩ(ξ;P ₁ ,P ₂), p=εE′Π(ξ;P ₁ ,P ₂), l=γ(P ₁ , P ₂)L  (30)With these definitions, we have introduced the scaled coordinateξ=x/l(t) (0≦ξ≦1), a small number ε(t), a length scale L(t) of the sameorder of magnitude as the fracture length l(t), and two dimensionlessevolution parameters P₁(t) and P₂(t), which depend monotonically on t.The form of the scaling (30) can be motivated from elementary elasticityconsiderations, by noting that the average aperture scaled by thefracture length is of the same order as the average net pressure scaledby the elastic modulus.

Four different scalings can be defined to emphasize above differentprimary limiting cases. These scalings yield power law dependence of L,ε, P₁, and P₂ on time t; i.e. L˜t^(α), ε˜t^(δ), P₁˜t^(β) ¹ , P₂˜t^(β) ², see Table 2 for the case of a radial fracture. Furthermore, theevolution parameters can take either the meaning of a toughness (K_(m),K_({tilde over (m)})), or a viscosity (M_(k), M _(k) ), or a storage(S_({tilde over (m)}), S _(k) ), or a leak-off coefficient (C_(m),C_(k)).

TABLE 2 Small parameter ε, lengthscale L, and parameters P₁ and P₂ forthe two storage scalings (viscosity and toughness) and the two leak-offscalings (viscosity and toughness). Scaling ε L P₁ P₂ storage/viscosity(M) $\left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}$$\left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/6}$$K_{m} = {K^{\prime}\left( \frac{1}{E^{\prime 3}\mu^{\prime}Q_{o}} \right)}^{1/4}$$C_{m} = {\text{C}^{\prime}\left( \frac{E^{\prime}t}{\mu^{\prime}Q_{o}^{3}} \right)^{1/6}}$storage/toughness (K)$\left( \frac{K^{\prime 4}}{E^{\prime 4}Q_{o}t} \right)^{1/3}$$\left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/3}$$M_{k} = {\mu^{\prime}\left( \frac{E^{\prime 3}Q_{o}}{K^{\prime 4}} \right)}$$C_{k} = {\text{C}^{\prime}\left( \frac{E^{\prime 4}t}{K^{\prime 4}Q_{o}^{2}} \right)^{1/6}}$leak-off/viscosity (M)$\left( \frac{\mu^{\prime}C^{\prime 2}}{E^{\prime}Q_{o}t} \right)^{1/4}$$\frac{Q_{o}t^{1/2}}{C^{\prime}}$ $K_{\overset{\sim}{m}} = K_{m}$$S_{\overset{\sim}{m}} = \left( \frac{\mu^{\prime}Q_{o}^{3}}{E^{\prime}C^{\prime 6}t} \right)^{1/4}$leak-off/toughness (K)$\left( \frac{K^{\prime 4}C^{\prime 2}}{E^{\prime 4}Q_{o}^{2}t} \right)^{1/4}$$\frac{Q_{o}t^{1/2}}{C^{\prime}}$ $M_{\overset{\sim}{k}} = M_{k}$$S_{\overset{\sim}{k}} = \left( \frac{K^{\prime 4}Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}t} \right)^{1/4}$

The regimes of solutions can be conceptualized in a rectangular phasediagram MK{tilde over (K)}{tilde over (M)} shown in FIG. 9. Each of thefour primary regimes (M, K, {tilde over (M)}, and {tilde over (K)}) ofhydraulic fracture propagation corresponding to the vertices of thediagram is dominated by only one component of fluid global balance whilethe other can be neglected (i.e. respective P₁=0, see Table 1) and onlyone dissipative process while the other can be neglected (i.e.respective P₂=0, see Table 1). As follows from the stationary tipsolution, the behavior of the solution at the tip also depends on theregime of solution: Ω˜(1−ρ)^(2/3) at the M-vertex, Ω˜(1−ρ)^(5/8) at the{tilde over (M)}-vertex, and Ω˜(1−ρ)^(1/2) at the K- and {tilde over(K)}-vertices.

The edges of the rectangular phase diagram MK{tilde over (K)}{tilde over(M)} can be identified with the four secondary limiting regimescorresponding to either the dominance of one of the two fluid globalbalance mechanisms or the dominance of one of the two energy dissipationmechanisms: storage-edge (MK, C_(m)=C_(k)=0), leak-off-edge ({tilde over(M)}{tilde over (K)}, S_({tilde over (m)})=S_({tilde over (k)})=0),viscosity-edge (M{tilde over (M)}, K_(m)=K_({tilde over (m)})=0), andK{tilde over (K)}-toughness-edge (M_(k)=M _(k) =0). The solution alongthe storage-edge MK and along the leak-off-edge {tilde over (M)}{tildeover (K)} has the property that it evolves with time t according to apower law, i.e., according to l˜t^(α) where the exponent α depends onthe regime of propagation: α=⅔ on the storage-edge MK and α=½ onleak-off-edge {tilde over (M)}{tilde over (K)}.

The regime of propagation evolves with time from the storage-edge to theleak-off edge since the parameters C's and S's depend on t, but not K'sand M's. With respect to the evolution of the solution in time, it isuseful to locate the position of the state point in the MK{tilde over(K)}{tilde over (M)} space in terms of η which is a power of any of theparameters K's and M's and a dimensionless time, eitherτ_(m{tilde over (m)})=t/t_(m{tilde over (m)}) or τ_(k k) =t/t_(k k)where

$\begin{matrix}{{\eta = \frac{K^{\prime 4}}{E^{\prime 3}\mu^{\prime}Q_{o}}},{t_{m\;\overset{\sim}{m}} = \frac{\mu^{\prime}Q_{o}^{3}}{E^{\prime}C^{\prime 6}}},{t_{k\overset{\sim}{k}} = \frac{K^{\prime 4}Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}}}} & (31)\end{matrix}$also noting that τ_(m{tilde over (m)})=ηt_(k{tilde over (k)}) since

$\begin{matrix}{\frac{t_{k\overset{\sim}{k}}}{t_{m\overset{\sim}{m}}} = \eta} & (32)\end{matrix}$The parameters M's, K's, C's and S's can be expressed in terms of η andτ_(m{tilde over (m)}) (or τ_(k k) ) according toK_(m)=K_({tilde over (m)})=η^(1/4), M_(k)=M_({tilde over (k)})=η⁻¹  (33)C_(m)=τ_(m{tilde over (m)}) ^(1/6), C_(k)=η^(−1/6)τ_(m{tilde over (m)})^(1/6)=τ_(k{tilde over (k)}) ^(1/6)  (34)S_({tilde over (m)})=τ_(m{tilde over (m)}) ^(1/4),S_({tilde over (k)})=η^(1/4)τ_(m{tilde over (m)})^(−1/4)=τ_(k{tilde over (k)}) ^(−1/4)  (35)

A point in the parametric space MK{tilde over (K)}{tilde over (M)} isthus completely defined by η and any of these two times. The evolutionof the state point can be conceptualized as moving along a trajectoryperpendicular to the storage- or the leak-off-edge.

In summary, the MK-edge corresponds to the origin of time, and the{tilde over (M)}{tilde over (K)}-edge to the end of time (except inimpermeable rocks). Thus, given all the problem parameters whichcompletely define the number η, the system evolves with time (e.g., timeτ_(mk)) along a η-trajectory, starting from the MK-edge (C_(m)=C_(k)=0)and ending at the {tilde over (M)}{tilde over (K)}-edge(S_({tilde over (k)})=S _(m) =0). If η=0, the fluid is inviscid (μ′=0)and the system then evolves along the toughness-edge from K to {tildeover (K)}. If η=∞, then K′=0 the system evolves along the viscosity-edgefrom M to {tilde over (M)}; The dependence of the scaled solution F canthus be expressed in the form F(ξ,τ;η), where τ is one of thedimensionless time, irrespective of the adopted scaling.

II. Embodiments Utilizing a Second Parametric Space

A. Radial Fractures

Determining the solution of the problem of a radial hydraulic fracturepropagating in a permeable rock consists of finding the aperture w ofthe fracture, and the net pressure p (the difference between the fluidpressure p_(f) and the far-field stress σ_(o)) as a function of both theradial coordinate r and time t, as well as the evolution of the fractureradius R(t). The functions R(t), w(r,t), and p(r,t) depend on theinjection rate Q_(o) and on the four material parameters E′, μ′, K′, andC′ respectively defined as

$\begin{matrix}{{E^{\prime} = \frac{E}{1 - v^{2}}}{\mu^{\prime} = {12\;\mu}}{K^{\prime} = {4\left( \frac{2}{\pi} \right)^{1/2}K_{Ic}}}{C^{\prime} = {2C_{I}}}} & (36)\end{matrix}$The three functions R(t), w(r,t), and p(r,t) are determined by solving aset of equations which can be summarized as follows.Elasticity Equation

$\begin{matrix}{w = {\frac{R}{E^{\prime}}{\int_{0}^{1}{{G\left( {{r/R},s} \right)}{p\left( {{sR},t} \right)}s\ {\mathbb{d}s}}}}} & (37)\end{matrix}$where G is a known elastic kernel. This singular integral equationexpresses the non-local dependence of the fracture width w on the netpressure p.Lubrication Equation

$\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{1}{r}\frac{\partial}{\partial r}\left( {{rw}^{3}\frac{\partial p}{\partial r}} \right)}} & (38)\end{matrix}$This non-linear differential equation governs the flow of viscousincompressible fluid inside the fracture. The function g(r,t) denotesthe rate of fluid leak-off, which evolves according to Carter's law

$\begin{matrix}{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(r)}}}} & (39)\end{matrix}$where t_(o)(r) is the exposure time of point r (i.e., the time at whichthe fracture front was at a distance r from the injection point).Global Volume Balance

$\begin{matrix}{{Q_{o}t} = {{2\pi{\int_{0}^{R}{w\ r{\mathbb{d}x}}}} + {2\pi{\int_{0}{r{\int_{0}^{R{(\tau)}}{{g\left( {r,\tau}\  \right)}{\mathbb{d}r}\ {\mathbb{d}\tau}}}}}}}} & (40)\end{matrix}$This equation expresses that the total volume of fluid pumped is equalto the sum of the fracture volume and the volume of fluid lost in therock surrounding the fracture.Propagation Criterion

$\begin{matrix}{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{R - r}}},{1 - {\frac{r}{R}{\operatorname{<<}1}}}} & (41)\end{matrix}$Within the framework of linear elastic fracture mechanics, this equationembodies fact that the fracture is always propagating and that energy isdissipated continuously in the creation of new surfaces in rock (at aconstant rate per unit surface)Tip Conditions

$\begin{matrix}{{w = 0},{{w^{3}\frac{\partial p}{\partial r}} = 0},{r = R}} & (42)\end{matrix}$The tip of the propagating fracture corresponds to a zero width and to azero fluid flow rate condition.1. Scalings

The general solution of this problem (which includes understanding thedependence of the solution on all the problem parameters) can beconsiderably simplified through the application of scaling laws. Scalingof this problem hinges on defining the dimensionless crack opening Ω,net pressure Π, and fracture radius γ asw=εLΩ(ρ;P ₁ ,P ₂), p=εE′Π(ρ;P ₁ ,P ₂), R=γ(P ₁ , P ₂)L  (43)

These definitions introduce the scaled coordinate ρ=r/R(t) (0≦ρ≦1), asmall number ε(t), a length scale L(t) of the same order of magnitude asthe fracture length R(t), and two dimensionless evolution parametersP₁(t) and P₂(t), which depend monotonically on t. As is shown below,three different scalings (“viscosity”, “toughness,” and “leak-off”) canbe defined, which yield power law dependence of L, ε, P₁, and P₂ on timet; i.e. L˜t^(α), ε˜t^(δ), P₁˜t^(β) ¹ , P₂˜t^(β) ² . The form of thescaling (43) can be motivated from elementary elasticity considerations,by noting that the average aperture scaled by the fracture radius is ofthe same order as the average net pressure scaled by the elasticmodulus.

The main equations are transformed as follows, under the scaling (43).

Elasticity Equation

$\begin{matrix}{\Omega = {\gamma{\int_{0}^{1}{{G\left( {\rho,s} \right)}{\Pi\left( {{s;P_{1}},P_{2}} \right)}s\ {\mathbb{d}s}}}}} & (44)\end{matrix}$Lubrication Equation

$\begin{matrix}{{{\left( {\frac{\overset{.}{ɛ}t}{ɛ} + \frac{\overset{.}{L}t}{L}} \right)\Omega} - {\frac{\overset{.}{L}t}{L}\rho\frac{\partial\Omega}{\partial\rho}} + {{\overset{.}{P}}_{1}^{t}\left( {\frac{\partial\Omega}{\partial P_{1}} - {\frac{\rho}{\gamma}\frac{\partial\gamma}{\partial P_{1}}\frac{\partial\Omega}{\partial\rho}}} \right)} + {{\overset{.}{P}}_{2}^{t}\left( {\frac{\partial\Omega}{\partial P_{2}} - {\frac{\rho}{\gamma}\frac{\partial\gamma}{\partial P_{2}}\frac{\partial\Omega}{\partial\rho}}} \right)} + {G_{c}\Gamma}} = {\frac{1}{G_{m}\gamma^{2}\rho}\frac{\partial}{\partial\rho}\left( {{\rho\Omega}^{3}\frac{\partial\Pi}{\partial\rho}} \right)}} & (45)\end{matrix}$where the leak-off function Γ(ρ; P₁, P₂) is defined as

${{\Gamma\left( {{\rho;P_{1}},P_{2}} \right)} = \frac{1}{\sqrt{1 - {t_{o}/t}}}},{t > t_{o}}$Global Mass Balance

$\begin{matrix}{{{2{\pi\gamma}^{2}{\int_{0}^{1}{{\Omega\rho}\ {\mathbb{d}\rho}}}} + {2{\pi G}_{c}{\int_{0}^{1}{u^{{2\alpha} - {1/2}}{\gamma^{2}\left( {{u^{\beta_{1}}P_{1}},{u^{\beta_{2}}P_{2}}} \right)}{I\left( {{u^{\beta_{1}}P_{1}},{u^{\beta_{2}}P_{2}}} \right)}\ {\mathbb{d}u}}}}} = G_{v}} & (46)\end{matrix}$where I is given by

I(X₁, X₂) = ∫₀¹Γ(ρ; X₁, X₂)ρ 𝕕ρPropagation CriterionΩ=G _(k)γ^(1/2)(1−ρ)^(1/2)1−ρ<<1  (47)Four dimensionless groups G_(ν), G_(m), G_(k), G_(c) appear in theseequations:

$\begin{matrix}{{G_{v} = \frac{Q_{o}t}{ɛ\; L^{3}}},{G_{m} = \frac{\mu^{\prime}}{ɛ^{3}E^{\prime}t}},{G_{k} = \frac{K^{\prime}}{ɛ\; E^{\prime}L^{1/2}}},{G_{c} = \frac{C^{\prime}t^{1/2}}{ɛ\; L}}} & (48)\end{matrix}$

While the group G_(ν) is associated with the volume of fluid pumped,G_(m), G_(k), and G_(c) can be interpreted as dimensionless viscosity,toughness, and leak-off coefficients, respectively. Three differentscalings can be identified, with each scaling leading to a differentdefinition of the set ε, L, P₁, and P₂. Each scaling is obtained bysetting G_(ν)=1 and one of the other groups to 1 (G_(m) for theviscosity scaling, G_(k) for the toughness scaling, and G_(c) for theleak-off scaling), with the two other groups being identified as P₁ andP₂. Three scalings denoted as viscosity, toughness, and leak-off canthus be defined depending on whether the group containing μ′ (G_(m)) K′(G_(k)) or C′ (G_(c)) is set to 1. The three scalings are summarized inTable 3.

TABLE 3 Small parameter ε, lengthscale L, and parameters P₁ and P₂ forthe viscosity, toughness, and leak-off scaling. Scaling ε L P₁ P₂Viscosity $\left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}$$\left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/9}$$K_{m} = {K^{\prime}\left( \frac{t^{2}}{\mu^{\prime 5}E^{\prime 13}Q_{o}^{3}} \right)}^{1/18}$$C_{m} = {\text{C}^{\prime}\left( \frac{E^{\prime 4}t^{7}}{\mu^{\prime 4}Q_{o}^{6}} \right)^{1/18}}$Toughness $\left( \frac{K^{\prime 6}}{E^{\prime 6}Q_{o}t} \right)^{1/5}$$\left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/5}$$M_{k} = {\mu^{\prime}\left( \frac{Q_{o}^{3}E^{\prime 13}}{K^{\prime 18}t^{2}} \right)}^{1/5}$$C_{k} = {\text{C}^{\prime}\left( \frac{E^{\prime 8}t^{3}}{K^{\prime 8}Q_{o}^{2}} \right)^{1/10}}$Leak-off $\left( \frac{C^{\prime 6}t}{Q_{o}^{2}} \right)^{1/4}$$\left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/4}$$K_{c} = {K^{\prime}\left( \frac{Q_{o}^{2}}{E^{\prime 8}C^{\prime 10}t^{3}} \right)}^{1/8}$$M_{c} = {\mu^{\prime}\mspace{11mu}\left( \frac{Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}t^{7}} \right)^{1/4}}$

The evolution of the radial fracture can be conceptualized in theternary phase diagram MKC shown in FIG. 10. First, however, thedimensionless number η and time τ are introduced as

$\begin{matrix}{{\eta = \frac{K^{\prime 14}}{E^{\prime 11}\mu^{'3}C^{\prime 4}Q_{o}}},{\tau = \frac{t}{t_{c\; m}}},{{{with}\mspace{14mu} t_{c\; m}} = \left( \frac{\mu^{\prime 4}Q_{o}^{6}}{E^{\prime 4}C^{\prime 18}} \right)^{1/7}}} & (49)\end{matrix}$

As shown in Table 3, the evolution parameters P₁ and P₂ in the threescalings can be expressed in terms of η and τ only. Both K_(m) and C_(m)are positive power of time τ, while K_(c) and M_(c) are negative powerof τ; furthermore, M_(k)˜τ^(−2/5) and C_(k)˜τ^(3/10). Hence, theviscosity scaling is appropriate for small time, while the leak-offscaling is appropriate for large time. The toughness scaling applies tointermediate time when both M_(κ) and C_(κ) are o(1).

TABLE 4 Dependence of the parameters P₁ and P₂ on the dimensionless timeτ and number η for the viscosity, toughness, and leak-off scaling.Scaling P₁ P₂ Viscosity K_(m) = η^(1/14)τ^(1/9) C_(m) = τ^(7/18)Toughness C_(k) = η^(−2/35)τ^(3/10) M_(k) = η^(−9/35)τ^(−2/5) Leak-offM_(c) = τ^(−7/4) K_(c) = η^(1/14)τ^(−3/8)The solution of a hydraulic fracture starts at the M-vertex (K_(m)=0,C_(m)=0) and ends at the C-vertex (M_(c)=0, K_(c)=0); it evolves withtime τ, along a trajectory which is controlled only by the number η, afunction of all the problem parameters (i.e., Q_(o), E′, μ′, K′, andC′). If η=0 (the rock has zero toughness), the evolution from M to C isdone directly along the base MC of the ternary diagram MKC. Withincreasing η (which can be interpreted for example as increasingrelative toughness, the trajectory is pulled towards the K-vertex. Forη=∞, two possibilities exist: either the rock is impermeable (C′=0) andthe system evolves directly from M to K, or the fluid is inviscid andthe system then evolves from K to C.

At each corner of the MKC diagram, there is only one dissipativemechanism at work; for example, at the M-vertex, energy is onlydissipated in viscous flow of the fracturing fluid since the rock isassumed to be impermeable and to have zero toughness. It is interestingto note that the mathematical solution is characterized by a differenttip singularity at each corner, reflecting the different nature of thedissipative mechanism.

M-cornerΩ˜(1−ρ)^(2/3)Π˜(1−ρ)^(−1/3) for ρ˜1  (50)C-cornerΩ˜(1−ρ)^(5/8)Π˜(1−ρ)^(−3/8) for ρ˜1  (51)K-cornerΩ˜(1−ρ)^(1/2)Π˜Const for ρ˜1  (52)

The transition of the solution in the tip region between two corners canbe analyzed by considering the stationary solution of a semi-infinitehydraulic fracture propagating at constant speed.

2. Applications of the Scaling Laws

The dependence of the scaled solution F={Ω, Π, γ} is thus of the formF(ρ,τ;η), irrespective of the adopted scaling. In other words, thescaled solution is a function of the dimensionless spatial and timecoordinates ρ and τ, which depends only on η, a constant for aparticular problem. Thus the laws of similitude between field andlaboratory experiments simply require that η is preserved and that therange of dimensionless time τ is the same—even for the general case whenneither the fluid viscosity, nor the rock toughness, nor the leak-off offracturing fluid in the reservoir can be neglected.

Although the solution in any scaling can readily be translated intoanother scaling, each scaling is useful because it is associated with aparticular process. Furthermore, the solution at a corner of the MKCdiagram in the corresponding scaling (i.e., viscosity at M, toughness atK, and leak-off at C) is self-similar. In other words, the scaledsolution at these vertices does not depend on time, which implies thatthe corresponding physical solution (width, pressure, fracture radius)evolves with time according to a power law. This property of thesolution at the corners of the MKC diagram is important, in part becausehydraulic fracturing near one corner is completely dominated by theassociated process. For example, in the neighborhood of the M-corner,the fracture propagates in the viscosity-dominated regime; in thisregime, the rock toughness and the leak-off coefficient can beneglected, and the solution in this regime is given for all practicalpurposes by the zero-toughness, zero-leak-off solution at the M-vertex.Findings from work along the MK edge where the rock is impermeablesuggest that the region where only one process is dominant is surprisinglarge. FIG. 11 shows the variation of γ_(mk) (the fracture radius in theviscosity scaling) with the dimensionless toughness K_(m) for animpermeable rock (K_(m)=0 corresponds to the M-vertex, K_(m)=∞ (i.e.,M_(c)=0) to the K-vertex). These results indicate that a hydraulicfracture propagating in an impermeable rock is in theviscosity-dominated regime if K_(m)<K_(mm)≅1, and in thetoughness-dominated regime if K_(m)>K_(mk)≅4.

Accurate solutions can be obtained for a radial hydraulic fracturepropagating in regimes corresponding to the edges MK, KC, and CM of theMKC diagram. These solutions enable one to identify the three regimes ofpropagation (viscosity, toughness, and leak-off).

The range of values of the evolution parameters P₁ and P₂ for which thefracture propagates in one of the primary regimes (viscosity, toughness,and leak-off) can be identified. The criteria in terms of the numbers P₁and P₂ can be translated in terms of the physical parameters (i.e., theinjection rate Q_(o), the fluid viscosity μ, the rock toughness K_(lc),the leak-off coefficient C_(l), and the rock elastic modulus E′).

The primary regimes of fracture propagation (corresponding to thevertices of the MKC diagram) are characterized by a simple power lawdependence of the solution on time. Along the edges of the MKC triangle,outside the regions of dominance of the corners, the evolution of thesolution can readily be tabulated.

In some embodiments of the present invention, the tabulated solutionsare used for quick design of hydraulic fracturing treatments. In otherembodiments, the tabulated solutions are used to interpret real-timemeasurements during fracturing, such as down-hole pressure.

The derived solutions can be considered as exact within the framework ofassumptions, since they can be evaluated to practically any desireddegree of accuracy. These solutions are therefore useful benchmarks totest numerical simulators currently under development.

3. Derivation of Solutions Along Edges of the Triangular ParametricSpace

Derivation of the solution along the edges of the triangle MKC and atthe C-vertex are now described. The identification of the differentregimes of fracture propagation are also described.

a. CK-Solution

Along the CK-edge of the MKC triangle, the influence of the viscosity isneglected and the solution depends only on one parameter (either K_(c),the dimensionless toughness in the leak-off scaling, or thedimensionless leak-off coefficient C_(k) in the toughness scalingC_(k)). In one embodiment, the solution is constructed starting from theimpermeable case (K-vertex) and it is evolved with increasing C_(k)towards the C-vertex.

Since the fluid is taken to be inviscid along the CK-edge, the pressuredistribution along the fracture is uniform and the corresponding openingis directly deduced by integration of the elasticity equation (44)

$\begin{matrix}{{\Pi_{kc} = {\Pi_{kc}\left( C_{k} \right)}},{\Omega_{kc} = {\frac{8}{\pi}\gamma_{kc}\Pi_{kc}\sqrt{1 - \rho^{2}}}}} & (53)\end{matrix}$Combining (53) with the propagation criterion (47) yields

$\begin{matrix}{{\Pi_{kc} = {\frac{\pi}{8\sqrt{2}}\gamma_{kc}^{{- 1}/2}}},{\Omega_{kc} = {\left( \frac{\gamma_{kc}}{2} \right)^{1/2}\sqrt{1 - \rho^{2}}}}} & (54)\end{matrix}$The radius γ_(kc) is determined as a function of C_(k). An equation forγ_(kc) can be deduced from the global balance of mass

$\begin{matrix}{{\frac{\mathbb{d}\gamma_{kc}}{\mathbb{d}C_{k}} = {{\frac{4}{3C_{k}}{\frac{\gamma_{k}^{5/2}}{\gamma_{kc}^{3/2}}\left\lbrack {1 - \left( \frac{\gamma_{kc}}{\gamma_{k}} \right)^{5/2}} \right\rbrack}} - {\frac{8\pi}{3}\frac{\gamma_{kc}^{1/2}}{\gamma_{k}^{5/2}}{I\left( C_{k} \right)}}}}{where}} & (55) \\{{{\gamma_{k} \equiv {\gamma_{kc}(0)}} = \left( \frac{3}{\pi\sqrt{2}} \right)^{5/2}},{{{and}\mspace{14mu}{I(X)}} = {\int_{0}^{1}{\frac{1}{\sqrt{1 - {\tau_{o}\left( {\rho,X} \right)}}}\rho\ {\mathbb{d}\rho}}}}} & (56)\end{matrix}$with τ_(o)=t_(o)(r)/t denoting the scaled exposure time of point r. Thefunction τ_(o)(ρ, X) can be found by inverting

$\begin{matrix}{\rho = {\tau_{o}^{2/5}\frac{\gamma_{kc}\left( {\tau_{o}^{3/10}X} \right)}{\gamma_{kc}(X)}}} & (57)\end{matrix}$which is deduced from the definition of ρ by taking into account thepower law dependence of L_(k) and C_(k) on time.

Since τ_(o) (1, X)=1, the integral I(X) defined in (56) is singular atρ=1. This singularity is weak, and its strength is known at X=0 and X=∞:X=0(τ_(o)=ρ^(5/2)) and at X=∞(τ_(o)=ρ⁴). From a computational point ofview, the integral can be calculated along the time axis with respect toτ_(o)

$\begin{matrix}\begin{matrix}{{I(X)} = {\frac{1}{\gamma_{kc}(X)}{\int_{0}^{1}{\frac{1}{{\tau_{o}^{3/5}\left( {1 - \tau_{o}} \right)}^{1/2}}\left\lbrack {{\frac{2}{5}\gamma_{kc}\left( {\tau_{o}^{3/10}X} \right)} +}\mspace{205mu} \right.}}}} \\{\left. {\frac{3}{10}\tau_{o}^{3/10}X\;{\gamma_{kc}^{\prime}\left( {\tau_{o}^{3/10}\ X} \right)}} \right\rbrack{\mathbb{d}\tau_{o}}}\end{matrix} & (58)\end{matrix}$

In some embodiments of the present invention, the solution can beobtained by solving the non-linear ordinary differential equation (55),using an implicit iterative algorithm.

b. MK-Solution

The MK-solution corresponds to regimes of fracture propagation inimpermeable rocks. One difficulty in obtaining this solution lies inhandling the changing nature of the tip behavior between the M- and theK-vertex. The tip asymptote is given by the classical square rootsingularity of linear elastic fracture mechanics (LEFM) wheneverK_(m)≠0. However, near the M-vertex (small K_(m)), the LEFM behavior isconfined to a small boundary layer, which does not influence thepropagation of the fracture. In this viscosity-dominated regime, thesingularity (50) develops as an intermediate asymptote.

The solution can be searched for in the form of a finite series of knownbase functions

$\begin{matrix}{\Pi_{k\; m} = {{\Pi_{o}^{*}\left( {\rho,M_{k}} \right)} + {\sum\limits_{i = 1}^{n_{\Pi}}{{A_{i}\left( M_{k} \right)}{\Pi_{i}^{*}(\rho)}}} + {{B\left( M_{k} \right)}{\Pi^{**}(\rho)}}}} & (59) \\{\;{{\overset{\_}{\Omega}}_{k\; m} = {{{\overset{\_}{\Omega}}_{o}^{*}\left( {\rho,M_{k}} \right)} + {\sum\limits_{i = 1}^{n_{\Omega}}{{C_{i}\left( M_{k} \right)}{{\overset{\_}{\Omega}}_{i}^{*}(\rho)}}} + {{B\left( M_{k} \right)}{{\overset{\_}{\Omega}}^{**}(\rho)}}}}} & (60)\end{matrix}$

where the introduction of Ω _(km)=Ω_(km)/γ_(km) excludes γ_(km) from theelasticity equation (44).

Since the non-linearity of the problem only arises from the lubricationequation (45), the series expansions (59) and (60) can be used tosatisfy the elasticity equation and the boundary conditions at the tipand at the inlet. In the proposed decomposition, the last terms {Π**,Ω**} are chosen such that the logarithmic pressure singularity near theinlet is satisfied. The corresponding opening is integrated bysubstituting this pressure function into (44). The first terms in theseries {Π_(o)*, Ω _(o)*} are constructed to exactly satisfy thepropagation equation and to account for the logarithmic pressureasymptote near the tip (which results from substituting the openingsquare root asymptote into the lubrication equation). It is alsorequired that {Π_(o)*, Ω _(o)*} exactly satisfy the elasticity equation(44). The regular part of the solution is represented by series of basefunctions {Π_(i)*, Ω _(i)*}. The choice of these functions is notunique; however, it seems consistent to require that Ω_(i)*˜(1−ρ)^(1/2+i) for ρ˜1. (The square root opening asymptote appearsonly in the first term, if one imposes that the function Π_(i)* does notcontribute to the stress intensity factor.) A convenient choice of thesebase functions are Jacobi polynomials with the appropriate weights.

Any pair {Π_(i)*, Ω _(i)*} does not satisfy the elasticity equation(44). Instead, the coefficients A_(i) and C_(i) are related by theelasticity equation through the matrix L_(ij) (which is independent ofM_(k) or K_(m)).

$\begin{matrix}{C_{i}^{({n_{\Omega},n_{\Pi}})} = {\sum\limits_{j = 1}^{n_{\Pi}}{L_{ij}A_{j}^{({n_{\Omega},n_{\Pi}})}}}} & (61)\end{matrix}$The problem is reduced to finding n_(Π)+1 unknown coefficients A_(i) andB, by solving the lubrication equation (45), which simplifies here to

$\begin{matrix}{{{{\rho\mspace{11mu}{\overset{\_}{\Omega}}_{km}^{3}\frac{\partial\Pi_{km}}{\partial\rho}} + {M_{k}{\int_{\rho}^{1}{{\overset{\_}{\Omega}}_{km}s\ {\mathbb{d}s}}}} + {\frac{4}{5}M_{k}\rho^{2}{\overset{\_}{\Omega}}_{km}} - {\frac{2}{5}{M_{k}^{2}\left\lbrack {{\int_{\rho}^{1}{\frac{\partial\;{\overset{\_}{\Omega}}_{km}}{\partial M_{k}}s\ {\mathbb{d}s}}} + {\frac{2}{\gamma_{km}}\frac{\mathbb{d}\gamma_{km}}{\mathbb{d}M_{k}}\left( {{\int_{\rho}^{1}\;{{\overset{\_}{\Omega}}_{km}s\ {\mathbb{d}s}}} + {\rho^{2}{\overset{\_}{\Omega}}_{km}}} \right)}} \right\rbrack}}} = 0}{{{where}\mspace{14mu}\gamma_{km}} = \left( {2\;\pi{\int_{\rho}^{1}\;{{\overset{\_}{\Omega}}_{km}\rho\ {\mathbb{d}\rho}}}} \right)^{{- 1}/3}}} & (62)\end{matrix}$

In some embodiments of the present invention, the lubrication equationis solved by an implicit iterative procedure. For example, the solutionat the current iteration can be found by a least squares method.

c. CM-Solution

In some embodiments, the solution along the CM-edge of the MKC triangleis found using the series expansion technique described above withreference to the MK-solution. In other embodiments, a numerical solutionis used based on the following algorithm.

The displacement discontinuity method is used to solve the elasticityequation (44). This method yields a linear system of equations betweenaperture and net pressure at nodes along the fracture. The coefficients(which can be evaluated analytically) need to be calculated only once asthey do not depend on C_(m). The lubrication equation (45) is solved bya finite difference scheme (either explicit or implicit). The fractureradius γ_(mc) is found from the global mass balance. Here, the numericaldifficulty is to calculate the amount of fluid lost due to the leak-off.

The propagation is governed by the asymptotic behavior of the solutionat the fracture tip. The tip asymptote can be used to establish arelationship between the opening at the computational node next to thetip and the tip velocity. However, this relationship evolves as C_(m)increases from 0 to ∞ (i.e., when moving from the M- to the C-vertex);it is obtained through a mapping of the autonomous solution of asemi-infinite hydraulic fracture propagating at constant speed in apermeable rock.

d. Solution Near the C-Vertex

The limit solution at the C-vertex, where both the viscosity and thetoughness are neglected, is degenerated as all the fluid injected intothe fracture has leaked into the rock. Thus the opening and the netpressure of the fracture is zero, while its radius is finite. In someembodiments of the present invention, the solution near the C-vertex isused for testing the numerical solutions along the CK and CM sides ofthe parametric triangle. The limitation of those solutions comes fromthe choice of the scaling. In order to approach the C-vertex, thecorresponding parameter (C_(k) or C_(m)) must grow indefinitely.Practically, these solutions are calculated up to some finite values ofthe parameters, for which they can be connected with asymptoticsolutions near the C-vertex along CM and CK sides. These asymptoticsolutions can be constructed as follows.

Consider first the CM-solution

F_(cm)={Ω_(cm)(ρ,M_(c)),Π_(cm)(ρ,M_(c)),γ_(cm)(M_(c))} near theC-vertex. It can be asymptotically approximated asγ_(cm)=γ_(c) +o(M _(c)), Ω_(cm) =M _(c) ^(α)γ_(c) Ω _(cm)(ρ)+o(M _(c)^(α)), Π_(cm) =M _(c) ^(α) Π _(cm)(ρ)+o(M _(c) ^(α)) (63)where γ_(c) denotes the finite fracture radius (in the leak-off scaling)at the C-vertex. The exponent α=¼ is determined by substituting theseexpansions into the lubrication equation (45), which then reduces to

$\begin{matrix}{\frac{\gamma_{c}}{\sqrt{1 - \rho^{4}}} = {\frac{1}{\rho}\frac{\mathbb{d}\;}{\mathbb{d}\rho}\left( {\rho\;{\overset{\_}{\Omega}}_{cm}^{3}\frac{\mathbb{d}{\overset{\_}{\Pi}}_{cm}}{\mathbb{d}\rho}} \right)}} & (64)\end{matrix}$

The asymptotic solution F _(cm)={ Ω _(cm)(ρ), Π _(cm)(ρ)} near theC-vertex is found by solving (64) along with the elasticity equation(44). This can be done using the series expansion technique describedabove. This problem is similar to the problem at the M-vertex (fracturepropagating in an impermeable solid with zero toughness), but with adifferent tip asymptote. Thus a set of base functions different from theone used for the M-corner are introduced.

The CK-solution F_(ck)={Ω_(ck)(ρ,K_(c)),Π_(ck)(ρ,K_(c)),γ_(ck)(K_(c))}near the C-vertex can also be sought in the form of an asymptoticexpansionγ_(ck)=γ_(c) +o(K _(c)), Ω_(ck) =K _(c) ^(β)γ_(c) Ω _(ck)(ρ)+o(K _(c)^(β)), Π_(ck) =K _(c) ^(β) Π _(ck)(ρ)+o(K _(c) ^(β))  (65)where β=1 is determined from the propagation condition (11). Thissolution is trivial, however, since the pressure is uniform; hence

$\begin{matrix}{\;{{{\overset{\_}{\Pi}}_{ck} = {\frac{\pi}{8}\left( {2\;\gamma_{c}} \right)^{{- 1}/2}}},{{\overset{\_}{\Omega}}_{ck} = {\left( {2\;\gamma_{c}} \right)^{{- 1}/2}\sqrt{1 - \rho^{2}}}}}} & (66)\end{matrix}$e. Regimes of Fracture Propagation

The regimes of fracture propagation can readily be identified once thesolutions at the vertices and along the edges of the MKC triangle havebeen tabulated using the algorithms and methods of solutions describedabove. Recall that for the parametric space under consideration, thereare three primary regimes of propagation (viscosity, toughness, andleak-off) associated with the vertices of the MKC triangle and that in acertain neighborhood of a corner, the corresponding process is dominant,see Table 5. For example, fracture propagation is in theviscosity-dominated regime if K_(m)<K_(mm) and C_(m)<C_(mm); in thisregion, the solution can be approximated for all practical purposes bythe zero-toughness, zero-leak-off solution at the M-corner (K_(m)=0,C_(m)=0).

TABLE 5 Range of the parameters P₁ and P₂ for which a primary process isdominant. Dominant Process Range on P₁ Range on P₂ Viscosity K_(m) <K_(mm) (M_(k) > M_(km)) C_(m) < C_(mm) (M_(c) > M_(cm)) Toughness C_(k)< C_(kk) (K_(c) > K_(ck)) M_(k) < M_(kk) (K_(m) > K_(mk)) Leak-off M_(c)< M_(cc) (C_(m) > C_(mc)) K_(c) < K_(cc) (C_(k) < C_(kc))

Identification of the threshold values of the evolution parameters (forexample, K_(mm) and C_(mm) for the viscosity-dominated regime) can beaccomplished by comparing the fracture radius with its reference valueat a corner. The corner process is considered as dominant, if thefracture radius is within 1% of its value at the corner. For example,K_(mm) and C_(mm) are deduced from the following conditions|γ_(mk)(K _(mm))−γ_(m)|/γ_(m)=1% |γ_(mk)(C _(mm))−γ_(m)|/γ_(m)=1%  (67)B. Plane Strain (KGD) Fractures1. Governing Equations and Boundary Conditions

$\begin{matrix}{Elasticity} & \; \\{{w = {\frac{l}{E^{\prime}}{\int_{0}^{1}{{G\left( {{x/l},s} \right)}{\Pi\left( {{sl},t} \right)}\ {\mathbb{d}s}}}}}{Lubrication}} & (68) \\{{\frac{\partial w}{\partial t} + g} = {\frac{1}{\mu^{\prime}}\frac{\partial\;}{\partial x}\left( {w^{3}\frac{\partial p}{\partial x}} \right)}} & (69)\end{matrix}$obtained by eliminating the radial flow rate q(x,t) between the fluidmass balance

$\begin{matrix}{{{\frac{\partial w}{\partial t} + \frac{\partial q}{\partial x} + g} = 0}{{and}\mspace{14mu}{Poiseuille}\mspace{14mu}{law}}} & (70) \\{{q = {{- \frac{w^{3}}{\mu^{\prime}}}\frac{\partial p}{\partial x}}}{{Leak}\text{-}{off}}} & (71) \\{g = \frac{C^{\prime}}{\sqrt{t - {t_{o}(x)}}}} & (72)\end{matrix}$where t_(o)(x) is the exposure time of point x

$\begin{matrix}{{Global}\mspace{14mu}{volume}\mspace{14mu}{balance}} & \; \\{{{Q_{o}t} = {{2{\int_{0}^{l{(t)}}{w\ {\mathbb{d}x}}}} + {2{\int_{0}^{t}{\int_{0}^{l{(\tau)}}{{g\left( {x,\tau} \right)}\ {\mathbb{d}x}\ {\mathbb{d}\tau}}}}}}}{{Propagation}\mspace{14mu}{criterion}}} & (73) \\{{{w \simeq {\frac{K^{\prime}}{E^{\prime}}\sqrt{l - x}}},{{1 - \frac{x}{l}} ⪡ 1}}{{Tip}\mspace{14mu}{conditions}}} & (74) \\{{w = 0},\mspace{14mu}{{w^{3}\frac{\partial p}{\partial x}} = 0},\mspace{14mu}{x = l}} & (75)\end{matrix}$2. Scaling

Similarly to the radial fracture, we define the dimensionless crackopening Ω, net pressure Π, and fracture length γ asw(x,t)=ε(t)L(t)Ω(ξ;P ₁ P ₂)  (76)p(x,t)=ε(t)E′Π(ξ;P ₁ ,P ₂)  (77)l(t)=γ(ξ;P ₁ P ₂)L(t)  (78)

These definitions introduce a scaled coordinate ξ=x/l(t) (0≦ξ≦1), asmall number ε(t), a length scale L(t) of the same order of magnitude asthe fracture length l(t), and two dimensionless parameters P₁(t), P₂(t)which depend monotonically on t. The form of the scaling (76)-(80) canbe motivated from elementary elasticity considerations, by noting thatthe average aperture scaled by the fracture radius is of the same orderas the average net pressure scaled by the elastic modulus. Explicitforms of the parameters ε(t), L(t), P₁(t), and P₂(t) are given below forthe viscosity, toughness, and leak-off scalings.

The main equations are transformed as follows, under the scaling(76)-(80).

Elasticity EquationΩ=γ∫₀ ¹ G(ξ,s)Π(s;P ₁ ,P ₂)ds  (79)Lubrication EquationThe left-hand side ∂w/∂t of the lubrication equation (69) can now bewritten as

$\begin{matrix}{{\frac{\partial w}{\partial t} + g} = {{\left( {{\overset{.}{ɛ}\; L} + {ɛ\;\overset{.}{L}}} \right)\Omega} - {ɛ\;\overset{.}{L}\;\xi\frac{\partial\Omega}{\partial\xi}} + {ɛ\; L{{\overset{.}{P}}_{1}\left( {\frac{\partial\Omega}{\partial P_{1}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{1}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {ɛ\; L{{\overset{.}{P}}_{2}\left( {\frac{\partial\Omega}{\partial P_{2}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{2}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {C^{\prime}t^{{- \; 1}/2}\Gamma\;\left( {{\xi;P_{1}},P_{2}} \right)}}} & (80)\end{matrix}$while the right hand side is transformed into

$\begin{matrix}{{\frac{1}{\mu^{\prime}}\frac{\partial\;}{\partial x}\left( {w^{3}\frac{\partial p}{\partial x}} \right)} = {\frac{ɛ^{4}E^{\prime}L}{\mu^{\prime}\gamma^{2}}\frac{\partial\;}{\partial\xi}\left( {\Omega^{3}\frac{\partial\Pi}{\partial\xi}} \right)}} & (81)\end{matrix}$The leak-off function Γ(ξ;P₁, P₂), which is defined as

$\begin{matrix}{{{\Gamma\;\left( {{\xi;P_{1}},P_{2}} \right)} = \frac{1}{1 - {t_{o}/t}}},{t > t_{o}}} & (82)\end{matrix}$can be computed as part of the solution, once the parameters P₁,P₂ havebeen identified. After multiplying both sides by t/εR, we obtain a newform of the lubrication equation

$\begin{matrix}{{{\left( {\frac{\overset{.}{ɛ}t}{ɛ} + \frac{\overset{.}{L}t}{L}} \right)\Omega} - {\frac{\overset{.}{L}t}{L}\xi\frac{\partial\Omega}{\partial\xi}} + {{\overset{.}{P}}_{1}{t\left( {\frac{\partial\Omega}{\partial P_{1}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{1}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {{\overset{.}{P}}_{2}{t\left( {\frac{\partial\Omega}{\partial P_{2}} - {\frac{\xi}{\gamma}\frac{\partial\gamma}{\partial P_{2}}\frac{\partial\Omega}{\partial\xi}}} \right)}} + {\frac{C^{\prime}t^{\;{1/2}}}{ɛ\; L}\Gamma}} = {\frac{ɛ^{3}E^{\prime}t}{\mu^{\prime}}\frac{1}{\gamma^{2}}\frac{\partial\;}{\partial\xi}\left( {\Omega^{3}\frac{\partial\Pi}{\partial\xi}} \right)}} & (83)\end{matrix}$Global Mass Balance

$\begin{matrix}{{{2\;\gamma{\int_{0}^{1}{\Omega\ {\mathbb{d}\rho}}}} + {2\frac{C^{\prime}t^{1/2}}{ɛ\; L}{\int_{0}^{1}{u^{a - {1/2}}{\gamma\left( {{\tau_{o}^{\beta_{1}}P_{1}},{\tau_{o}^{\beta_{2}}P_{2}}} \right)}{I\left( {{\tau_{o}^{\beta_{1}}P_{1}},{\tau_{o}^{\beta_{2}}P_{2}}} \right)}\ {\mathbb{d}u}}}}} = \frac{Q_{o}t}{ɛ\; L^{2}}} & (84)\end{matrix}$where I is given by I(X₁,X₂)=∫₀ ¹Γ(ξ;X₁,X₂)dξPropagation Criterion

$\begin{matrix}{\Omega = {{{\frac{K^{\prime}}{ɛ\; E^{\prime}L^{1/2}}{\gamma^{1/2}\left( {1 - \xi} \right)}^{1/2}\mspace{14mu} 1} - \xi} ⪡ 1}} & (85)\end{matrix}$These equations show that there are 4 dimensionless groups: G_(ν),C_(m), C_(k), G_(c) (only G_(ν) differs from the radial case, in view ofthe different dimension of Q_(o))

$\begin{matrix}{{G_{v} = \frac{Q_{o}t}{ɛ\; L^{2}}},{G_{m} = \frac{\mu^{\prime}}{ɛ^{3}E^{\prime}t}},{G_{k} = \frac{K^{\prime}}{ɛ\; E^{\prime}L^{1/2}}},{G_{c} = \frac{C^{\prime}t^{1/2}}{ɛ\; L}}} & (86)\end{matrix}$a. Viscosity Scaling.

The small parameter ε_(m) and the lengthscale L_(m) are determined bysetting G_(ν)=1 and G_(m)=1. Hence,

$\begin{matrix}{{ɛ_{m} = \left( \frac{\mu^{\prime}}{E^{\prime}t} \right)^{1/3}},{L_{m} = \left( \frac{E^{\prime}Q_{o}^{3}t^{4}}{\mu^{\prime}} \right)^{1/6}}} & (87)\end{matrix}$The two parameters P₁=G_(k) and P₂=G₂ are identified as K_(m) and C_(m),a dimensionless toughness and a dimensionless leak-off coefficient,respectively

$\begin{matrix}{{K_{m} = {K^{\prime}\left( \frac{1}{E^{\prime 3}\mu^{\prime}Q_{o}} \right)}^{1/4}},{C_{m} = {C^{\prime}\left( \frac{E^{\prime}t}{\mu^{\prime}Q_{o}^{3}} \right)}^{1/6}}} & (88)\end{matrix}$b. Toughness Scaling.Now, ε_(k) and L_(k) are determined from G_(ν)=1 and G_(k)=1. Hence,

$\begin{matrix}{{ɛ_{k} = \left( \frac{K^{\prime 4}}{E^{\prime 4}Q_{o}t} \right)^{1/2}},{L_{k} = \left( \frac{E^{\prime}Q_{o}t}{K^{\prime}} \right)^{2/3}}} & (89)\end{matrix}$The two parameters P₁=G_(m) and P₂=G_(c) correspond to M_(k) and C_(k),a dimensionless viscosity and a dimensionless leak-off coefficient,respectively

$\begin{matrix}{{M_{k} = {\mu^{\prime}\left( \frac{E^{\prime 3}Q_{o}}{K^{\prime 4}} \right)}},{C_{k} = {C^{\prime}\left( \frac{E^{\prime 4}t}{K^{\prime 4}Q_{o}^{2}} \right)}^{1/6}}} & (90)\end{matrix}$c. Leak-off Scaling.

Finally, the leak-off scaling corresponds to G_(ν)=1 and G_(c)=1. Hence,

$\begin{matrix}{ɛ_{c} = {{\frac{C^{\prime 2}}{Q_{o}}\mspace{14mu} L_{c}} = \left( \frac{Q_{o}^{2}t}{C^{\prime 2}} \right)^{1/2}}} & (91)\end{matrix}$and the two parameters P₁=G_(k) and P₂=G_(m) are now identified as K_(c)and M_(c), a dimensionless viscosity and a dimensionless toughness,respectively

$\begin{matrix}{{K_{c} = {K^{\prime}\left( \frac{Q_{o}^{2}}{E^{\prime 4}C^{\prime 6}t} \right)}^{1/4}},{M_{c} = {\mu^{\prime}\left( \frac{Q_{o}^{3}}{E^{\prime}C^{\prime 6}t} \right)}}} & (92)\end{matrix}$

We note that both C_(k), C_(m) are positive power of time t while K_(c),M_(c) are negative power of t. Hence, the leak-off scaling isappropriate for large time, and either the viscosity scaling or thetoughness scaling is appropriate for small time. As discussed below, thesolution starts from a point on the MK-side of a ternary parameter space(C_(k)=0, C_(m)=0) and tends asymptotically towards the C-point(M_(c)=0, K_(c)=0), following a straight trajectory which is controlledby a certain number η, a function of all the problem parameters exceptC′ (i.e., Q_(o), E′, μ′, K′).

3. Time Scales

It is of interest to express the small parameters ε's, the length scalesL's, and the dimensionless parameters M's, K's, and C's in terms of timescales. Two time scales t_(m), t_(k), are naturally defined as

$\begin{matrix}{{t_{m} = \frac{\mu^{\prime}}{E^{\prime}}},{t_{k} = \frac{K^{\prime 4}}{E^{\prime 4}Q_{o}}}} & (93)\end{matrix}$Note that unlike the radial fracture, it is not possible to define acharacteristic time t_(c), since Q_(o) has the dimension squared of C′.Hence,

$\begin{matrix}{{ɛ_{m} = \left( \frac{t_{m}}{t} \right)^{1/3}},{ɛ_{k} = \left( \frac{t_{k}}{t} \right)^{1/3}}} & (94) \\{{L_{m} = {\left( \frac{t}{t_{m}} \right)^{2/3}{\overset{\_}{L}}_{m}}},{L_{k} = {\left( \frac{t}{t_{k}} \right)^{2/3}{\overset{\_}{L}}_{k}}}} & (95)\end{matrix}$where the L's are intrinsic length scales defined as

$\begin{matrix}{{{\overset{\_}{L}}_{m} = \left( \frac{\mu^{\prime\; 3}Q_{o}^{3}}{E^{\prime 3}} \right)^{1/6}},{{\overset{\_}{L}}_{k} = \left( \frac{K^{\prime}}{E^{\prime}} \right)^{2}}} & (96)\end{matrix}$

Next, consider the dimensionless parameters M's, K's, and C's which canbe rewritten in terms of the characteristic times t_(cm), and t_(ck)

$\begin{matrix}{{C_{m} = \left( \frac{t}{t_{cm}} \right)^{1/6}},{C_{k} = \left( \frac{t}{t_{ck}} \right)}} & (97) \\{{{M_{c} = \left( \frac{t_{cm}}{t} \right)},{K_{c} = \left( \frac{t_{ck}}{t} \right)^{1/4}}}{where}} & (98) \\{{t_{cm} = {{ɛ_{c}^{- 3}t_{m}} = \left( \frac{\mu^{\prime}Q_{o}^{3}}{E^{\prime}C^{\prime 6}} \right)}},{t_{ck} = {{ɛ_{c}^{- 3}t_{k}} = \left( \frac{K^{\prime\; 4}Q_{o}^{2}}{E^{\prime\; 4}C^{\prime 6}} \right)}}} & (99)\end{matrix}$

It is thus convenient to introduce a parameter η related to the ratiosof characteristic times, which is defined as

$\begin{matrix}{\eta = \frac{K^{\prime 4}}{E^{\prime 3}\mu^{\prime}Q_{o}}} & (100)\end{matrix}$

Indeed, it is easy to show that the various characteristic time ratioscan be expressed in terms of η

$\begin{matrix}{\frac{t_{ck}}{t_{cm}} = \eta} & (101)\end{matrix}$

Note also that η can be expressed as

$\begin{matrix}{\eta = \frac{{\overset{\_}{L}}_{k}^{2}}{{\overset{\_}{L}}_{m}^{2}}} & (102)\end{matrix}$

Furthermore, if we introduce the dimensionless time τ

$\begin{matrix}{\tau = \frac{t}{t_{cm}}} & (103)\end{matrix}$(acknowledging at the same time that the choice of t_(cm) to scale thetime is arbitrary, as t_(ck) could have been used as well), theparameters M's, K's, and C's can be expressed in terms of τ and η asfollowsK_(m)=η^(1/4), C_(m)=τ^(1/6), C_(k)=η^(−1/6)τ^(1/6)  (104)M_(k)=η⁻¹, M_(c)=τ⁻¹, K_(c)=η^(1/4)τ^(−1/4)  (105)

The dependence of the scaled solution F={Ω,Π,γ} is thus of the formF(ξ,τ;η), irrespective of the adopted scaling (but γ=γ(τ;η)). In otherwords, the scaled solution is a function of dimensionless spatial andtime coordinate, ξ and τ, which depends on only one parameter, η, whichis constant for a particular problem. Thus the laws of similitudebetween field and laboratory experiments simply require that η ispreserved and that the range of dimensionless time τ is the same—evenfor the general case of viscosity, toughness, and leak-off.

It is again convenient to introduce the ternary diagram MKC shown inFIG. 12. With time τ, the system evolves along a η-trajectory (alongwhich η is a constant), starting from a point on the MK-side and endingat the C-vertex. If η=0 (the rock has zero toughness), the evolutionfrom M to C is done directly along the base BC of the ternary diagramMKC. For η=∞, the fluid is inviscid and the system then evolves from Kto C.

The KGD fracture differs from the radial fracture by the existence ofonly characteristic time rather than two for the penny-shaped fracture.The characteristic number η for the KGD fracture is independent of theleak-off coefficient C′, which only enters the scaling of time.

4. Relationship Between Scalings

Any scaling can be translated into any of the other two. It can readilybe established that

$\begin{matrix}{{{K_{m} = M_{k}^{{- 1}/4}},{C_{m} = M_{c}^{{- 1}/6}},{C_{k} = K_{c}^{{- 2}/3}}}{and}} & (106) \\{\frac{ɛ_{m}}{ɛ_{k}} = {M_{k}^{1/3} = K_{m}^{{- 4}/3}}} & (107) \\{\frac{ɛ_{k}}{ɛ_{c}} = {K_{c}^{4/3} = C_{k}^{- 2}}} & (108) \\{\frac{ɛ_{c}}{ɛ_{m}} = {C_{m}^{2} = M_{c}^{{- 1}/3}}} & (109) \\{\frac{L_{m}}{L_{k}} = {M_{k}^{{- 1}/6} = K_{m}^{2/3}}} & (110) \\{\frac{L_{k}}{L_{c}} = {K_{c}^{{- 2}/3} = C_{k}}} & (111) \\{\frac{L_{c}}{L_{m}} = {C_{m}^{- 1} = M_{c}^{1/6}}} & (112)\end{matrix}$III. Applications

Applications of hydraulic fracturing include the recovery of oil and gasfrom underground reservoirs, underground disposal of liquid toxic waste,determination of in-situ stresses in rock, and creation of geothermalenergy reservoirs. The design of hydraulic fracturing treatmentsbenefits from information that characterize the fracturing fluid, thereservoir rock, and the in-situ state of stress. Some of theseparameters are easily determined (such as the fluid viscosity), but forothers, it is more difficult (such as physical parameters characterizingthe reservoir rock and in-situ state of stress).

By utilizing the various embodiments of the present invention, the“difficult” parameters can be assessed from measurements (such asdownhole pressure) collected during a hydraulic fracturing treatment.The various embodiments of the present invention recognize that scaledmathematical solutions of hydraulic fractures with simple geometrydepend on only two numbers that lump time and all the physicalparameters describing the problem. There are many different ways tocharacterize the dependence of the solution on two numbers, as describedin the different sections above, and all of these are within the scopeof the present invention.

Various parametric spaces have been described, and trajectories withinthose spaces have also been described. Each trajectory shows a pathwithin the corresponding parametric space that describes the evolutionof a particular treatment over time for a given set of physicalparameter values. That is to say, each trajectory lumps all of thephysical parameters, except time. Since there exists a unique solutionat each point in a given parametric space, which needs to be calculatedonly once and which can be tabulated, the evolution of the fracture canbe computed very quickly using these pre-tabulated solutions. In someembodiments, pre-tabulated points are very close together in theparametric space, and the closest pre-tabulated point is chosen as asolution. In other embodiments, solutions are interpolated betweenpre-tabulated points.

The various parametric spaces described above are useful to performparameter identification, also referred to as “data inversion.” Datainversion involves solving the so-called “forward model” many times,where the forward model is the tool to predict the evolution of thefracture, given all the problems parameters. Data inversion alsoinvolves comparing predictions from the forward model with measurements,to determine the set of parameters that provide the best match betweenpredicted and measured responses.

Historically, running forward models has been computationally demanding,especially given the complexity of the hydraulic fracturing process.Utilizing the various embodiments of the present invention, however, theforward model includes pre-tabulated scaled solutions in terms of twodimensionless parameters, which only need to be “unscaled” throughtrivial arithmetic operations. These developments, and others, makepossible real-time, or near real-time, data inversion while performing ahydraulic fracturing treatment.

Although the present invention has been described in conjunction withcertain embodiments, it is to be understood that modifications andvariations may be resorted to without departing from the spirit andscope of the invention as those skilled in the art readily understand.Such modifications and variations are considered to be within the scopeof the invention and the appended claims. For example, the scope of theinvention encompasses the so-called power law fluids (a generalizedviscous fluid characterized by two parameters and which degenerates intoa Newtonian fluid when the power law index is equal to 1). Also forexample, the scope of the invention encompasses the evolution of thehydraulic fracture following “shut-in” (when the injection of fluid isstopped). Hence, various embodiments of the invention contemplateinterpreting data gathered after shut-in.

1. A method of hydraulic fracturing comprising: performing a singleinjection of a viscous fluid wherein injecting the viscous fluidcomprises injecting the viscous fluid at a substantially constantvolumetric rate; measuring a pressure of the viscous fluid; determininga value of at least one dimensionless parameter associated with thepressure; and determining a value of a physical parameter from the atleast one dimensionless parameter to monitor the fracture.
 2. The methodof claim 1 wherein at least one of the dimensionless parametersrepresents a toughness of reservoir rock.
 3. The method of claim 1wherein at least one of the dimensionless parameters represents a fluidstorage mechanism.
 4. The method of claim 1 wherein at least one of thedimensionless parameters represents a dimensionless leak-offcoefficient.
 5. The method of claim 1, wherein the at least onedimensionless parameter represents a dimensionless leak-off coefficient.6. The method of claim 1, wherein the physical parameter comprises afracture dimension.
 7. The method of claim 1, wherein the at least onedimensionless parameter comprises a dimensionless crack opening.
 8. Themethod of claim 1, wherein the at least one dimensionless parametercomprises a dimensionless fracture radius.